senior-capstone.Rmd

---
title: "Senior capstone project"
output: html_notebook
---

This is the executable code for my senior capstone project, "Investigating the roles of innate immune cells in intestinal immunity and gut microbiome signaling." 

---
Pre-processing
---

First, the relevant packages are loaded. 

```{r Loading packages}

library(Seurat)
library(patchwork)
library(pathfindR)
library(sqldf)
library(dplyr)
library(ggplot2)
library(cowplot)
library(radiant.data)
library(org.Hs.eg.db)
library(textshape)
library(KEGGgraph)

```

The relevant files from indicated datasets (GSE150050, GSE185224, and GSE125527) are imported from GEO into the R environment. They are then converted to expression matrices. 

GSE150050 contains scRNA-seq data of CD127+ cells (ILCs) from four tissue locations. This data is narrowed using SQL on metadata files to isolate colon samples and clustered ('Seurat clustering & ILC phenotype analysis'); ILC expression data is be extracted and subsequently used in the 'Metabolic pathway analysis' section. 

```{r Importing datasets (GSE 150050)}

## Read in the counts data
counts.15 <- read.csv("C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/Data/GSE150050/STAR_raw_counts.csv", row.names=1)

## Read in the metadata
metadata.15 <- read.csv("C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/Data/GSE150050/GSE150050_metadata.csv", sep = ";", row.names=1)

## SQL script is used to narrow metadata.15 to colon samples
metadata15 <- as.data.frame(metadata.15)
metadata.15.colon <- sqldf("SELECT * FROM metadata15 WHERE TISSUE = 'COLON'")

## Metadata indicates all colon samples have distinct identifier 'GNI' in cell ID
counts15 <- as.data.frame(counts.15)
counts.15.colon <- counts15[,grepl("GNI", colnames(counts15))]

## Seurat object creation
seurat.15 <- CreateSeuratObject(counts = counts.15)

```

GSE185224 contains primary scRNA-seq data from the small intestines and colon epithelium of three donors. This data is be clustered ('Seurat clustering & ILC phenotype analysis'); ILC expression data is extracted and subsequently used in the 'Metabolic pathway analysis' section.

```{r Importing datasets (GSE 185224)}

## Read in the counts data
counts.18 <- Read10X_h5("C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/Data/GSE185224/GSE185224_Donor1_filtered_feature_bc_matrix.h5", use.names=TRUE)

## Seurat object creation from Gene Expression subset (not Antibody Capture subset)
seurat.18 <- CreateSeuratObject(counts = counts.18$`Gene Expression`)

```

GSE125527 contains primary scRNA-seq data of immune cells taken from fifteen patients presenting with or without ulcerative colitis/IBD. This data is divided by condition, processed, and re-integrated; it is then clustered with condition factored into analysis ('Condition-based clustering and phenotype composition analysis'); ILC expression data is extracted and subsequently used in the 'Metabolic pathway analysis' section.

```{r Importing datasets (GSE 125527)}

## Intestinal immune cells counts data from directory and append
dir.12.int <- "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/Data/GSE125527/UMI/HealthyI/"
files.12.int <- list.files(path = dir.12.int, pattern = ".tsv.gz", full.names = TRUE)

### Read in and aggregate  samples
counts.12.int <- read.csv(files.12.int[1],sep="\t", row.names=1)
for(i in 2:length(files.12.int)){
  counts.12b.int <- read.csv(files.12.int[i],sep="\t", row.names=1)
  counts.12.int <- rbind(counts.12.int, counts.12b.int)
}
counts.12.int <- t(counts.12.int)

## Seurat objects creation
seurat.12.int <- CreateSeuratObject(counts = counts.12.int)


## PBMC counts data from directory and append
dir.12.healthy <- "C:/Users/maiabennett/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/Data/GSE125527/UMI/HealthyPBMC/"
dir.12.UC <- "C:/Users/maiabennett/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/Data/GSE125527/UMI/UCPBMC/"
files.12.pbmc.healthy <- list.files(path = dir.12.healthy, pattern = ".tsv.gz", full.names = TRUE)
files.12.pbmc.UC <- list.files(path = dir.12.UC, pattern = ".tsv.gz", full.names = TRUE)

### Read in and aggregate healthy control samples
counts.12.healthy <- read.csv(files.12.pbmc.healthy[1],sep="\t", row.names=1)
for(i in 2:length(files.12.pbmc.healthy)){
  counts.12b.healthy <- read.csv(files.12.pbmc.healthy[i],sep="\t", row.names=1)
  counts.12.healthy <- rbind(counts.12.healthy, counts.12b.healthy)
}
counts.12.pbmc.healthy <- t(counts.12.healthy)

### Make metadata for later analysis
metadata.12.healthy <- data.frame(x1= colnames(counts.12.pbmc.healthy), x2 = "healthy")
colnames(metadata.12.healthy) <- c("barcode", "condition")

### Read in and aggregate UC samples
counts.12.UC <- read.csv(files.12.pbmc.UC[1],sep="\t", row.names=1)
for(i in 2:length(files.12.pbmc.UC)){
  counts.12b.UC <- read.csv(files.12.pbmc.UC[i],sep="\t", row.names=1)
  counts.12.UC <- rbind(counts.12.UC, counts.12b.UC)
}
counts.12.pbmc.UC <- t(counts.12.UC)

### Make metadata for later analysis
metadata.12.UC <- data.frame(x1= colnames(counts.12.pbmc.UC), x2 = "UC")
colnames(metadata.12.UC) <- c("barcode", "condition")
metadata.12.UC <- column_to_rownames(metadata.12.UC, loc = 1)

## Seurat objects creation
seurat.12.healthy <- CreateSeuratObject(counts = counts.12.pbmc.healthy, metadata = metadata.12.healthy)
seurat.12.UC <- CreateSeuratObject(counts = counts.12.pbmc.UC, metadata = metadata.12.UC)

## Create combined object
counts.12 <- rbind(t(counts.12.pbmc.healthy), t(counts.12.pbmc.UC))
counts.12 <- t(counts.12)
metadata.12 <- rbind(metadata.12.healthy, metadata.12.UC)
metadata.12 <- column_to_rownames(metadata.12, loc = 1)
seurat.12 <- CreateSeuratObject(counts = counts.12)
seurat.12 <- AddMetaData(seurat.12, metadata.12, col.name = "condition")

```

Then, the data pre-processed for further analysis. 

```{r Converting datasets to Seurat objects (GSE 150050)}

## QC and valid cell selection
seurat.15[["percent.mt"]] <- PercentageFeatureSet(seurat.15, pattern = "^MT-")
VlnPlot(seurat.15, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
seurat.15 <- subset(seurat.15, subset = nFeature_RNA > 200 & nFeature_RNA < 6000)

## Normalization of data
seurat.15 <- NormalizeData(seurat.15)

### Feature selection
seurat.15 <- FindVariableFeatures(seurat.15, selection.method = "vst", nfeatures = 2000)

### Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(seurat.15), 10)

### Plot variable features with and without labels
plot1 <- VariableFeaturePlot(seurat.15)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat15-varfeatures.png", width= 900, height=480)
plot1
plot2
dev.off()

## Scaling of data
all.genes <- rownames(seurat.15)
seurat.15 <- ScaleData(seurat.15, features = all.genes)

## Perform linear dimension reduction
seurat.15 <- RunPCA(seurat.15, features = VariableFeatures(object = seurat.15))
print(seurat.15[["pca"]], dims = 1:5, nfeatures = 5)
DimPlot(seurat.15, reduction = "pca")
DimHeatmap(seurat.15, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(seurat.15, dims = 1:15, cells = 500, balanced = TRUE)
seurat.15 <- JackStraw(seurat.15, num.replicate = 100)
seurat.15 <- ScoreJackStraw(seurat.15, dims = 1:20)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat15-jackstraw.png", width= 700, height=480)
JackStrawPlot(seurat.15, dims = 1:15)
dev.off()

## Save object for future recall
saveRDS(seurat.15, file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat15-preprocessed.rds")

```


```{r Converting datasets to Seurat objects (GSE 185224)}

## QC and valid cell selection
seurat.18[["percent.mt"]] <- PercentageFeatureSet(seurat.18, pattern = "^MT-")
VlnPlot(seurat.18, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
seurat.18 <- subset(seurat.18, subset = nFeature_RNA > 200 & nFeature_RNA < 6000)

## Normalization of data
seurat.18 <- NormalizeData(seurat.18)

### Feature selection
seurat.18 <- FindVariableFeatures(seurat.18, selection.method = "vst", nfeatures = 2000)

### Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(seurat.18), 10)

### Plot variable features with and without labels
plot3 <- VariableFeaturePlot(seurat.18)
plot4 <- LabelPoints(plot = plot3, points = top10, repel = TRUE)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat18-varfeatures.png", width= 900, height=480)
plot3 + plot4
dev.off()

## Scaling of data
all.genes <- rownames(seurat.18)
seurat.18 <- ScaleData(seurat.18, features = all.genes)

## Perform linear dimension reduction
seurat.18 <- RunPCA(seurat.18, features = VariableFeatures(object = seurat.18))
print(seurat.18[["pca"]], dims = 1:5, nfeatures = 5)
DimPlot(seurat.18, reduction = "pca")
DimHeatmap(seurat.18, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(seurat.18, dims = 1:15, cells = 500, balanced = TRUE)
seurat.18 <- JackStraw(seurat.18, num.replicate = 100)
seurat.18 <- ScoreJackStraw(seurat.18, dims = 1:20)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat18-jackstraw.png", width= 700, height=480)
JackStrawPlot(seurat.18, dims = 1:15)
dev.off()

## Save object for future recall
saveRDS(seurat.18, file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat18-preprocessed.rds")

```


```{r Converting datasets to Seurat objects (GSE 125527)}

# Intestinal immune cells
## QC and valid cell selection
seurat.12.int[["percent.mt"]] <- PercentageFeatureSet(seurat.12.int, pattern = "^MT-")
VlnPlot(seurat.12.int, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
seurat.12.int <- subset(seurat.12.int, subset = nFeature_RNA > 200 & nFeature_RNA < 6000)

## Normalization of data
seurat.12.int <- NormalizeData(seurat.12.int)

### Feature selection
seurat.12.int <- FindVariableFeatures(seurat.12.int, selection.method = "vst", nfeatures = 2000)

### Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(seurat.12.int), 10)

### Plot variable features with and without labels
plot1 <- VariableFeaturePlot(seurat.12.int)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12int-varfeatures.png", width= 900, height=480)
plot1 + plot2
dev.off()

## Scaling of data
all.genes <- rownames(seurat.12.int)
seurat.12.int <- ScaleData(seurat.12.int, features = all.genes)

## Perform linear dimension reduction
seurat.12.int <- RunPCA(seurat.12.int, features = VariableFeatures(object = seurat.12.int))
print(seurat.12.int[["pca"]], dims = 1:5, nfeatures = 5)
DimPlot(seurat.12.int, reduction = "pca")
DimHeatmap(seurat.12.int, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(seurat.12.int, dims = 1:15, cells = 500, balanced = TRUE)
seurat.12.int <- JackStraw(seurat.12.int, num.replicate = 100)
seurat.12.int <- ScoreJackStraw(seurat.12.int, dims = 1:20)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12int-jackstraw.png", width= 700, height=480)
JackStrawPlot(seurat.12.int, dims = 1:15)
dev.off()

## Save object for future recall
saveRDS(seurat.12.int, file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat12int-preprocessed.rds")


# PBMCs by condition
seurat.12.healthy[["percent.mt"]] <- PercentageFeatureSet(seurat.12.healthy, pattern = "^MT-")
VlnPlot(seurat.12.healthy, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
seurat.12.healthy <- subset(seurat.12.healthy, subset = nFeature_RNA > 200 & nFeature_RNA < 4000)

seurat.12.UC[["percent.mt"]] <- PercentageFeatureSet(seurat.12.UC, pattern = "^MT-")
VlnPlot(seurat.12.UC, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
seurat.12.UC <- subset(seurat.12.UC, subset = nFeature_RNA > 200 & nFeature_RNA < 3000)

## Normalization of data
seurat.12.healthy <- NormalizeData(seurat.12.healthy)
seurat.12.UC <- NormalizeData(seurat.12.UC)

### Feature selection
seurat.12.healthy <- FindVariableFeatures(seurat.12.healthy, selection.method = "vst", nfeatures = 2000)
seurat.12.UC <- FindVariableFeatures(seurat.12.UC, selection.method = "vst", nfeatures = 2000)

top10 <- head(VariableFeatures(seurat.12.healthy), 10)
top10 <- head(VariableFeatures(seurat.12.UC), 10)

### Plot individual variable features with and without labels
print(plot5 <- VariableFeaturePlot(seurat.12.healthy))
print(plot6 <- LabelPoints(plot = plot5, points = top10, repel = TRUE))
print(plot7 <- VariableFeaturePlot(seurat.12.UC))
print(plot8 <- LabelPoints(plot = plot7, points = top10, repel = TRUE))
png(file = "C:/Users/maiabennett/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/12healthy-varfeatures.png", width= 900, height=480)
plot5 + plot6
dev.off()
png(file = "C:/Users/maiabennett/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12UC-varfeatures.png", width= 900, height=480)
plot7 + plot8
dev.off()

### Normalize and identify variable features in combined Seurat object
seurat.12.list <- SplitObject(seurat.12, split.by = "condition")
seurat.12.list <- lapply(X = seurat.12.list, FUN = function(x) {
  x <- NormalizeData(x)
  x <- FindVariableFeatures(x, selection.method = "vst", nfeatures=2000)
})

### Select features that are repeatedly variable across datasets for integration
seurat.12.features <- SelectIntegrationFeatures(object.list = seurat.12.list)
seurat.12.anchors <- FindIntegrationAnchors(object.list = seurat.12.list, anchor.features = seurat.12.features)

### Create an 'integrated' data assay
seurat.12 <- IntegrateData(anchorset = seurat.12.anchors)
DefaultAssay(seurat.12) <- "integrated"

### Scaling and linear dimension reduction
seurat.12 <- ScaleData(seurat.12, verbose = FALSE)
seurat.12 <- RunPCA(seurat.12, npcs = 30, verbose = FALSE)
print(seurat.12[["pca"]], dims = 1:5, nfeatures = 5)
DimHeatmap(seurat.12, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(seurat.12, dims = 1:15, cells = 500, balanced = TRUE)
seurat.12 <- JackStraw(seurat.12, num.replicate = 100)
seurat.12 <- ScoreJackStraw(seurat.12, dims = 1:20)
png(file = "C:/Users/maiabennett/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12-jackstraw.png", width= 700, height=480)
JackStrawPlot(seurat.12, dims = 1:15)
dev.off()

seurat.12.healthy <- ScaleData(seurat.12.healthy, verbose = FALSE)
seurat.12.healthy <- RunPCA(seurat.12.healthy, npcs = 30, verbose = FALSE)
print(seurat.12.healthy[["pca"]], dims = 1:5, nfeatures = 5)
                                                                                                                                                                                                      
seurat.12.UC <- ScaleData(seurat.12.UC, verbose = FALSE)
seurat.12.UC <- RunPCA(seurat.12.UC, npcs = 30, verbose = FALSE)
print(seurat.12.UC[["pca"]], dims = 1:5, nfeatures = 5)

## Save object for future recall
saveRDS(seurat.12.healthy, file = "C:/Users/maiabennett/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat12healthy-preprocessed.rds")
## Save object for future recall
saveRDS(seurat.12.UC, file = "C:/Users/maiabennett/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat12UC-preprocessed.rds")
## Save object for future recall
saveRDS(seurat.12, file = "C:/Users/maiabennett/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat12-preprocessed.rds")

```

---
Clustering and phenotype composition analysis
---

This section contains code for Seurat clustering & ILC phenotype assignment.

```{r Initial clustering and cell type assignment (GSE 150050)}

## Cluster the cells
seurat.15 <- FindNeighbors(seurat.15, dims = 1:10)
seurat.15 <- FindClusters(seurat.15, resolution = 0.5)

### Look at cluster IDs of the first 5 cells
head(Idents(seurat.15), 5)

## Run non-linear dimensional reduction (UMAP/tSNE)
seurat.15 <- RunUMAP(seurat.15, dims = 1:10)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat15-initialUMAP.png")
DimPlot(seurat.15, reduction = "umap", label = TRUE)
dev.off()

## Find markers that differentiate clusters from one another (to help assign identities)
seurat.15.markers <- FindAllMarkers(seurat.15, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
manual.curate.15 <- seurat.15.markers %>%
  group_by(cluster) %>%
  slice_max(n = 3, order_by = avg_log2FC)
write.table(manual.curate.15, file="manual.curate.15.txt", row.names=TRUE, col.names=TRUE)

## Visualize specific markers to identify immune cells 
### ILC1s: CD127/IL7R+, CD161/KLRB1+, CD117/KIT-, CD294/PTGDR2-, CD3-
### ILC2s: CD127/IL7R+, CD161/KLRB1+, CD294/PTGDR2+, GATA3+, CD3-
### ILC3s: CD127/IL7R+, CD161/KLRB1+, CD117/KIT+, CD3-
### NKs: CD56/NCAM1+,CD3-, EOMES+
### Ts: CD3/CD3D+ (definite)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat15-initialmarkers.png")
FeaturePlot(seurat.15, features = c("CD3D", "NCAM1", "EOMES", "IL7R", "PTGDR2", "KIT", "KLRB1", "GATA3"), min.cutoff = "q9")
dev.off()

## Assign cell type identity to clusters
new.cluster.ids <- c("ILC3", "T cell", "ILC1", "ILC3", "Indeterminate", "Indeterminate", "Indeterminate ILC", "T cell", "T cell", "NK cell", "Indeterminate", "ILC2", "Indeterminate", "Indeterminate ILC")
names(new.cluster.ids) <- levels(seurat.15)
seurat.15 <- RenameIdents(seurat.15, new.cluster.ids)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat15-labeledUMAP.png")
DimPlot(seurat.15, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
dev.off()

## Visualize specific markers by cluster
#Idents(seurat.15) <- factor(Idents(seurat.15), levels = c("ILC3", "T cell", "ILC1", "ILC3", "Indeterminate", "Indeterminate", "Indeterminate ILC", "T cell", "T cell", "NK cell", "Indeterminate", "ILC2", "Indeterminate", "Indeterminate ILC"))
markers.to.plot <- c("CD3D", "NCAM1", "EOMES", "IL7R", "PTGDR2", "KIT", "KLRB1", "GATA3")
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat15-dotplot.png", width= 700, height=480)
DotPlot(seurat.15, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8) +
  RotatedAxis()
dev.off()

## Subset to cells with determinate and non-T cell identities
seurat.15.ilc <- subset(seurat.15, idents = c("ILC1", "ILC2", "ILC3", "NK cell", "Indeterminate ILC"))

## Cluster subset
seurat.15.ilc <- FindNeighbors(seurat.15.ilc, dims = 1:10)
seurat.15.ilc <- FindClusters(seurat.15.ilc, resolution = 0.5)

## Run non-linear dimensional reduction (UMAP/tSNE)
seurat.15.ilc <- RunUMAP(seurat.15.ilc, dims = 1:10)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat15ilc-initialUMAP.png")
DimPlot(seurat.15.ilc, reduction = "umap", label = TRUE)
dev.off()

## Find markers that differentiate clusters from one another (to help assign identities)
seurat.15.ilc.markers <- FindAllMarkers(seurat.15.ilc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
manual.curate.15.ilc <- seurat.15.ilc.markers %>%
  group_by(cluster) %>%
  slice_max(n = 3, order_by = avg_log2FC)
write.table(manual.curate.15.ilc, file="manual.curate.15.ilc.txt", row.names=TRUE, col.names=TRUE)

## Visualize specific markers to identify immune cell types
### ILC1s: CD127/IL7R+, CD161/KLRB1+, CD117/KIT-, CD294/PTGDR2-
### ILC2s: CD127/IL7R+, CD161/KLRB1+, CD117/KIT+-, CD294/PTGDR2+
### ILC3s: CD127/IL7R+, CD161/KLRB1+, CD117/KIT+, CD294/PTGDR-
### NKs: CD56/NCAM1+, EOMES+
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat15ilc-initialmarkers.png", width= 700, height=480)
FeaturePlot(seurat.15.ilc, features = c("NCAM1", "EOMES", "IL7R", "PTGDR2", "KIT", "KLRB1", "GATA3"), min.cutoff = "q9")
dev.off()

## Assign cell type identity to clusters
new.cluster.ids <- c("ILC3", "ILC3", "ILC1", "ILC3", "NK", "NK", "ILC2", "ILC2", "Indeterminate")
names(new.cluster.ids) <- levels(seurat.15.ilc)
seurat.15.ilc <- RenameIdents(seurat.15.ilc, new.cluster.ids)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat15ilc-labeledUMAP.png")
DimPlot(seurat.15.ilc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
dev.off()

## Visualize specific markers by cell type
#Idents(seurat.15.ilc) <- factor(Idents(seurat.15.ilc), levels = c("ILC3", "ILC3", "ILC1", "ILC3", "NK", "NK", "ILC2", "ILC2", "Indeterminate"))
markers.to.plot <- c("NCAM1", "EOMES", "IL7R", "PTGDR2", "KIT", "KLRB1", "GATA3")
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat15ilc-dotplot.png")
DotPlot(seurat.15.ilc, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8) +
  RotatedAxis()
dev.off()

## Save objects for future recall
saveRDS(seurat.15, file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat15-clustered.rds")
saveRDS(seurat.15.ilc, file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat15ilc-clustered.rds")

```


```{r Initial clustering and cell type assignment (GSE 185224)}

## Cluster the cells
seurat.18 <- FindNeighbors(seurat.18, dims = 1:10)
seurat.18 <- FindClusters(seurat.18, resolution = 0.5)

### Look at cluster IDs of the first 5 cells
head(Idents(seurat.18), 5)

## Run non-linear dimensional reduction (UMAP/tSNE)
seurat.18 <- RunUMAP(seurat.18, dims = 1:10)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat18-initialUMAP.png")
DimPlot(seurat.18, reduction = "umap", label = TRUE)
dev.off()

## Find markers that differentiate clusters from one another (to help assign identities)
seurat.18.markers <- FindAllMarkers(seurat.18, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
manual.curate.18 <- seurat.18.markers %>%
  group_by(cluster) %>%
  slice_max(n = 2, order_by = avg_log2FC)
write.table(manual.curate.18, file="manual.curate.18.txt", row.names=TRUE, col.names=TRUE)

## Assign cell type identity to clusters
new.cluster.ids <- c("Distal enterocyte", "Indeterminate enterocyte", "Undifferentiated", "Proximal enterocyte", "Proximal enterocyte", "Paneth", "Undifferentiated", "Intestinal globlet", "Undifferentiated", "Distal enterocyte", "Immune", "Proximal enterocyte", "Indeterminate enterocyte", "Paneth", "Immune", "Paneth", "Enteroendocrine")
names(new.cluster.ids) <- levels(seurat.18)
seurat.18 <- RenameIdents(seurat.18, new.cluster.ids)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat18-labeledUMAP.png")
DimPlot(seurat.18, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
dev.off()

## Subset clusters with immune identity
seurat.18.immune <- subset(seurat.18, idents = "Immune")

## Cluster immune cell subset
seurat.18.immune <- FindNeighbors(seurat.18.immune, dims = 1:10)
seurat.18.immune <- FindClusters(seurat.18.immune, resolution = 0.5)

### Look at cluster IDs of the first 5 cells
head(Idents(seurat.18.immune), 5)

## Run non-linear dimensional reduction (UMAP/tSNE)
seurat.18.immune <- RunUMAP(seurat.18.immune, dims = 1:10)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat18immune-initialUMAP.png")
DimPlot(seurat.18.immune, reduction = "umap", label = TRUE)
dev.off()

## Find markers that differentiate clusters from one another (to help assign identities)
seurat.18.immune.markers <- FindAllMarkers(seurat.18.immune, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
manual.curate.18.immune <- seurat.18.immune.markers %>%
  group_by(cluster) %>%
  slice_max(n = 2, order_by = avg_log2FC)
write.table(manual.curate.18.immune, file="manual.curate.18.immune.txt", row.names=TRUE, col.names=TRUE)

## Visualize specific markers to identify immune cells from a predominant epithelial population
## ILCs (CD127/IL7R+, CD117/KIT +, CD161/KLRB1+); no distinct signatures found
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat18immune-markers.png", width= 700, height=480)
FeaturePlot(seurat.18.immune, features = c("CD4", "CD8A", "IL7R", "PTGDR2", "KIT", "CD14", "CD19", "KLRB1"), min.cutoff = "q9")
dev.off()

# new.cluster.ids <- c("")
# names(new.cluster.ids) <- levels(seurat.18.immune)
# seurat.18 <- RenameIdents(seurat.18, new.cluster.ids)
# DimPlot(seurat.18.immune, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
# 
# ## Plot specific markers
# Idents(seurat.18.immune) <- factor(Idents(seurat.18.immune), levels = c(""))
# markers.to.plot <- c("CD4", "CD8A", "IL7R", "PTGDR2", "KIT", "CD14", "CD19", "KLRB1")
# DotPlot(seurat.18.immune, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8) +
#   RotatedAxis()

# Save objects for future recall
saveRDS(seurat.18, file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat18-clustered.rds")
saveRDS(seurat.18.immune, file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat18immune-clustered.rds")

```

This section contains code for the intestinal immune cell samples from GSE125527. This code is modeled on other samples from the same dataset; however, it has not been tested yet.

```{r Initial clustering and cell type assignment (GSE 125527)}

# Intestinal immune cells
## Cluster the cells
seurat.12.int <- FindNeighbors(seurat.12.int, dims = 1:10)
seurat.12.int <- FindClusters(seurat.12.int, resolution = 0.5)

### Look at cluster IDs of the first 5 cells
head(Idents(seurat.12.int), 5)

## Run non-linear dimensional reduction (UMAP/tSNE)
seurat.12.int <- RunUMAP(seurat.12.int, dims = 1:10)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12int-initialUMAP.png")
DimPlot(seurat.12.int, reduction = "umap", label = TRUE)
dev.off()

## Find markers that differentiate clusters from one another (to help assign identities)
seurat.12.int.markers <- FindAllMarkers(seurat.12.int, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
manual.curate.12.int <- seurat.12.int.markers %>%
  group_by(cluster) %>%
  slice_max(n = 3, order_by = avg_log2FC)
write.table(manual.curate.12.int, file="manual.curate.12.int.txt", row.names=TRUE, col.names=TRUE)

## Visualize specific markers to identify immune cells 
### ILC1s: CD127/IL7R+, CD161/KLRB1+, CD117/KIT-, CD294/PTGDR2-, CD3-, ID2+, CD8-, TBX21+
### ILC2s: CD127/IL7R+, CD161/KLRB1+, CD294/PTGDR2+, GATA3+, CD3-, ID2+, CD8-
### ILC3s: CD127/IL7R+, CD161/KLRB1+, CD117/KIT+, CD3-, ID2+, CD8-, IL23R+
### NKs: CD56/NCAM1+,CD3-, EOMES+, TBX21+
### Ts: CD3/CD3D+ (definite)
### Used CellMarker2.0 to call myeloids and B cells
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12int-initialmarkers.png", width= 700, height=480)
FeaturePlot(seurat.12.int, features = c("CD3D", "CD8A", "CD14", "CD19", "NCAM1", "EOMES", "IL7R",  "KLRB1", "GATA3", "ID2", "IL23R", "TBX21"), min.cutoff = "q9")
dev.off()

## Assign cell type identity to clusters
new.cluster.ids <- c("Dendritic cell", "Indeterminate", "Indeterminate", "Indeterminate", "B cell", "NK or T cell", "B cell", "Dendritic cell", "NK or T cell", "B cell", "Indeterminate", "Dendritic cell", "Dendritic cell")
names(new.cluster.ids) <- levels(seurat.12.int)
seurat.12.int <- RenameIdents(seurat.12.int, new.cluster.ids)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12int-labeledUMAP.png")
DimPlot(seurat.12.int, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
dev.off()

## Visualize specific markers by cluster
markers.to.plot <- c("CD3D", "CD8A", "CD14", "CD19", "NCAM1", "EOMES", "IL7R",  "KLRB1", "GATA3", "ID2", "IL23R", "TBX21")
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12int-dotplot.png")
DotPlot(seurat.12.int, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8) +
  RotatedAxis()
dev.off()

## Subset to cells with possible ILC marker expression and non-determinate T cell, B cell, or dendritic cell identity
seurat.12.int.inter <- subset(seurat.12.int, idents = c("Indeterminate", "NK or T cell"))

## Cluster subset
seurat.12.int.inter <- FindNeighbors(seurat.12.int.inter, dims = 1:10)
seurat.12.int.inter <- FindClusters(seurat.12.int.inter, resolution = 0.5)

## Run non-linear dimensional reduction (UMAP/tSNE)
seurat.12.int.inter <- RunUMAP(seurat.12.int.inter, dims = 1:10)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12intinter-initialUMAP.png")
DimPlot(seurat.12.int.inter, reduction = "umap", label = TRUE)
dev.off()

## Find markers that differentiate clusters from one another (to help assign identities)
seurat.12.int.inter.markers <- FindAllMarkers(seurat.12.int.inter, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
manual.curate.12.int.inter <- seurat.12.int.inter.markers %>%
  group_by(cluster) %>%
  slice_max(n = 3, order_by = avg_log2FC)
write.table(manual.curate.12.int.inter, file="manual.curate.12.int.inter.txt", row.names=TRUE, col.names=TRUE)

## Visualize specific markers to identify immune cells 
### ILC1s: CD127/IL7R+, CD161/KLRB1+, CD117/KIT-, CD294/PTGDR2-, CD3-, ID2+, CD8-, TBX21+
### ILC2s: CD127/IL7R+, CD161/KLRB1+, CD294/PTGDR2+, GATA3+, CD3-, ID2+, CD8-
### ILC3s: CD127/IL7R+, CD161/KLRB1+, CD117/KIT+, CD3-, ID2+, CD8-, IL23R+
### NKs: CD56/NCAM1+,CD3-, EOMES+, TBX21+
### Ts: CD3/CD3D+ (definite)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12intinter-initialmarkers.png", width= 700, height=480)
FeaturePlot(seurat.12.int.inter, features = c("CD3D", "CD8A", "CD14", "CD19", "NCAM1", "EOMES", "IL7R",  "KLRB1", "GATA3", "ID2", "IL23R", "TBX21"), min.cutoff = "q9")
dev.off()

## Subset to cells not expressing non-ILC marker CD8
seurat.12.int.inter2 <- subset(seurat.12.int.inter, idents = c(0,1,2,6,7))

## Cluster subset
seurat.12.int.inter2 <- FindNeighbors(seurat.12.int.inter2, dims = 1:10)
seurat.12.int.inter2 <- FindClusters(seurat.12.int.inter2, resolution = 0.5)

## Run non-linear dimensional reduction (UMAP/tSNE)
seurat.12.int.inter2 <- RunUMAP(seurat.12.int.inter2, dims = 1:10)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12intinter2-initialUMAP.png")
DimPlot(seurat.12.int.inter2, reduction = "umap", label = TRUE)
dev.off()

## Find markers that differentiate clusters from one another (to help assign identities)
seurat.12.int.inter2.markers <- FindAllMarkers(seurat.12.int.inter2, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
manual.curate.12.int.inter2 <- seurat.12.int.inter2.markers %>%
  group_by(cluster) %>%
  slice_max(n = 3, order_by = avg_log2FC)
write.table(manual.curate.12.int.inter2, file="manual.curate.12.int.inter2.txt", row.names=TRUE, col.names=TRUE)

## Visualize specific markers to identify immune cells 
### ILC1s: CD127/IL7R+, CD161/KLRB1+, CD117/KIT-, CD294/PTGDR2-, CD3-, ID2+, CD8-, TBX21+
### ILC2s: CD127/IL7R+, CD161/KLRB1+, CD294/PTGDR2+, GATA3+, CD3-, ID2+, CD8-
### ILC3s: CD127/IL7R+, CD161/KLRB1+, CD117/KIT+, CD3-, ID2+, CD8-, IL23R+
### NKs: CD56/NCAM1+,CD3-, EOMES+, TBX21+
### Ts: CD3/CD3D+ (definite)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12intinter2-initialmarkers.png", width= 700, height=480)
FeaturePlot(seurat.12.int.inter, features = c("CD3D", "CD8A", "CD14", "CD19", "NCAM1", "EOMES", "IL7R",  "KLRB1", "GATA3", "ID2", "IL23R", "TBX21"), min.cutoff = "q9")
dev.off()

## Subset to cells with possible ILC identity based on CD127 and CD161 expression
seurat.12.int.ilc <- subset(seurat.12.int.inter, idents = c(0,1,3,4,5))

## Cluster subset
seurat.12.int.ilc <- FindNeighbors(seurat.12.int.ilc, dims = 1:10)
seurat.12.int.ilc <- FindClusters(seurat.12.int.ilc, resolution = 0.5)

## Run non-linear dimensional reduction (UMAP/tSNE)
seurat.12.int.ilc <- RunUMAP(seurat.12.int.ilc, dims = 1:10)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12intilc-initialUMAP.png")
DimPlot(seurat.12.int.ilc, reduction = "umap", label = TRUE)
dev.off()

## Find markers that differentiate clusters from one another (to help assign identities)
seurat.12.int.ilc.markers <- FindAllMarkers(seurat.12.int.ilc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
manual.curate.12.int.ilc <- seurat.12.int.ilc.markers %>%
  group_by(cluster) %>%
  slice_max(n = 3, order_by = avg_log2FC)
write.table(manual.curate.12.int.ilc, file="manual.curate.12.int.ilc.txt", row.names=TRUE, col.names=TRUE)

## Visualize specific markers to identify immune cells 
### ILC1s: CD127/IL7R+, CD161/KLRB1+, CD117/KIT-, CD294/PTGDR2-, CD3-, ID2+, CD8-, TBX21+
### ILC2s: CD127/IL7R+, CD161/KLRB1+, CD294/PTGDR2+, GATA3+, CD3-, ID2+, CD8-
### ILC3s: CD127/IL7R+, CD161/KLRB1+, CD117/KIT+, CD3-, ID2+, CD8-, IL23R+
### NKs: CD56/NCAM1+,CD3-, EOMES+, TBX21+
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12intilc-initialmarkers.png", width= 700, height=480)
FeaturePlot(seurat.12.int.ilc, features = c("CD3D", "CD8A", "CD14", "CD19", "NCAM1", "EOMES", "IL7R",  "KLRB1", "GATA3", "ID2", "IL23R", "TBX21"), min.cutoff = "q9")
dev.off()

## Assign cell type identity to clusters
## CD3A was prevalent across all clusters, even after multiple runs of subsetting
## And ILC determining markers are not found in the matrix (can also be expressed on T cells)
## As such, ILCs will be assesed in groups likely containing ILCs based on marker expression (1,2)
## This is likely due to the 'mirroring' of Th subsets by ILCs
new.cluster.ids <- c("Non-ILCs", "Contains ILCs", "Contains ILCs", "Non-ILCs", "Non-ILCs")
names(new.cluster.ids) <- levels(seurat.12.int.ilc)
seurat.12.int.ilc <- RenameIdents(seurat.12.int.ilc, new.cluster.ids)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12intilc-labeledUMAP.png")
DimPlot(seurat.12.int.ilc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
dev.off()

## Visualize specific markers by cell type
#Idents(seurat.12.int.ilc) <- factor(Idents(seurat.12.int.ilc), levels = c(""))
markers.to.plot <- c("CD3D", "CD8A", "CD14", "CD19", "NCAM1", "EOMES", "IL7R",  "KLRB1", "GATA3", "ID2", "IL23R", "TBX21")
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12intilc-dotplot.png")
DotPlot(seurat.12.int.ilc, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8) +
  RotatedAxis()
dev.off()

## Save objects for future recall
saveRDS(seurat.12.int, file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat12int-clustered.rds")
saveRDS(seurat.12.int.inter, file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat12intinter-clustered.rds")
saveRDS(seurat.12.int.inter2, file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat12intinter2-clustered.rds")
saveRDS(seurat.12.int.ilc, file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat12intilc-clustered.rds")

```

This section contains code for composition analysis of expression in datasets (graphics, comparisons, etc).

```{r Comparing phenotype composition}

## Visualize specific microbiome interaction genes by condition and cell type
## GSE 150050
markers.to.plot <- c("NOD2", "CARD9", "ATG16L1*", "IRGM")
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat15-microdotplot.png")
DotPlot(seurat.15.ilc, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8) +
  RotatedAxis()
dev.off()

seurat.15.ilc$celltype <- Idents(seurat.15.ilc)
plots <- VlnPlot(seurat.15.ilc, features = c("NOD2", "CARD9", "ATG16L1", "IRGM"), group.by = "celltype",
    pt.size = 0, combine = FALSE, y.max = 2)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat15ilc-vinmicrogenes.png")
wrap_plots(plots = plots, ncol = 1)
dev.off()


## GSE 125527 intestines
markers.to.plot <- c("NOD2", "CARD9", "ATG16L1*", "IRGM")
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12intilc-microdotplot.png")
DotPlot(seurat.12.int.ilc, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8) +
  RotatedAxis()
dev.off()

seurat.12.int.ilc$celltype <- Idents(seurat.12.int.ilc)
plots <- VlnPlot(seurat.12.int.ilc, features = c("NOD2", "CARD9", "ATG16L1", "IRGM"), group.by = "celltype",
    pt.size = 0, combine = FALSE, y.max = 2)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12intilc-vinmicrogenes.png")
wrap_plots(plots = plots, ncol = 1)
dev.off()

```

---
Condition-based clustering and phenotype composition analysis
---

This section contains code for the integrated analysis of GSE125527 by condition (healthy or UC). 

```{r Integrated clustering and cell type assignment (GSE 125527)}

## Cluster the cells
seurat.12 <- FindNeighbors(seurat.12, dims = 1:10)
seurat.12 <- FindClusters(seurat.12, resolution = 0.5)

### Look at cluster IDs of the first 5 cells
head(Idents(seurat.12), 5)

## Run non-linear dimensional reduction (UMAP/tSNE)
seurat.12 <- RunUMAP(seurat.12, dims = 1:10)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12-initialUMAP.png")
DimPlot(seurat.12, reduction = "umap", label = TRUE)
dev.off()
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12-initialUMAPcond.png")
DimPlot(seurat.12, reduction = "umap", group.by = "condition")
dev.off()

## Find markers that differentiate clusters from one another (to help assign identities)
DefaultAssay(seurat.12) <- "RNA"
seurat.12.markers <- FindAllMarkers(seurat.12, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
manual.curate.12 <- seurat.12.markers %>%
  group_by(cluster) %>%
  slice_max(n = 2, order_by = avg_log2FC)
write.table(manual.curate.12, file="manual.curate.12.txt", row.names=TRUE, col.names=TRUE)

## Visualize specific markers to identify cell types
### ILC1s: CD127/IL7R+, CD161/KLRB1+, ID2+, T-bet/TBX21+, CD294/PTGDR2-, CD117/KIT-
## ILC2s: CD127/IL7R+, CD161/KLRB1+, ID2+, CD294/PTGDR2+, GATA3+, CD117/KIT+-
### ILC3s: CD127/IL7R+, CD161/KLRB1+, ID2+, CD294/PTGDR-, CD117/KIT+, IL23R+
### NKs: CD56/NCAM1+, EOMES+
### Ts: CD3/CD3D+ (definite), TBX21+
### Monocytes: CD14+
### Neutrophils, Basophils, Eosinophils, DCs: CD11b+
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12-initialmarkers.png", width= 700, height=480)
FeaturePlot(seurat.12, features = c("CD3D", "CD8A", "ITGAM", "CD14", "CD19", "NCAM1", "EOMES", "IL7R", "KLRB1", "ID2", "IL2RA", "THY1", "TBX21", "GATA3", "IL23R"), min.cutoff = "q9")
dev.off()


## Subset to cells with possible ILC identity (large bunch of clusters with scattered ILC marker expression, general lymphoid identity markers)
seurat.12.inter <- subset(seurat.12, idents = c(0, 1, 2, 4, 6, 7, 8, 9, 14))

## Cluster the cells
DefaultAssay(seurat.12.inter) <- "integrated"
seurat.12.inter <- FindNeighbors(seurat.12.inter, dims = 1:10)
seurat.12.inter <- FindClusters(seurat.12.inter, resolution = 0.5)

### Look at cluster IDs of the first 5 cells
head(Idents(seurat.12.inter), 5)

## Run non-linear dimensional reduction (UMAP/tSNE)
seurat.12.inter <- RunUMAP(seurat.12.inter, dims = 1:10)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12inter-initialUMAP.png")
DimPlot(seurat.12.inter, reduction = "umap", label = TRUE)
dev.off()
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12inter-initialUMAPcond.png")
DimPlot(seurat.12.inter, reduction = "umap", group.by = "condition")
dev.off()

## Find markers that differentiate clusters from one another (to help assign identities)
DefaultAssay(seurat.12.inter) <- "RNA"
seurat.12.inter.markers <- FindAllMarkers(seurat.12.inter, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
manual.curate.12.ilc <- seurat.12.inter.markers %>%
  group_by(cluster) %>%
  slice_max(n = 2, order_by = avg_log2FC)
write.table(manual.curate.12.ilc, file="manual.curate.12.ilc.txt", row.names=TRUE, col.names=TRUE)

## Visualize specific markers to identify cell types
### ILC1s: CD127/IL7R+, CD161/KLRB1+, ID2+, T-bet/TBX21+, CD294/PTGDR2-, CD117/KIT-
## ILC2s: CD127/IL7R+, CD161/KLRB1+, ID2+, CD294/PTGDR2+, GATA3+, CD117/KIT+-
### ILC3s: CD127/IL7R+, CD161/KLRB1+, ID2+, CD294/PTGDR-, CD117/KIT+, IL23R+
### NKs: CD56/NCAM1+, EOMES+
### Ts: CD3/CD3D+ (definite), TBX21+
### Monocytes: CD14+
### Neutrophils, Basophils, Eosinophils, DCs: CD11b+
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12inter-initialmarkers.png", width= 700, height=480)
FeaturePlot(seurat.12.inter, features = c("CD3D", "CD8A", "ITGAM", "CD14", "CD19", "NCAM1", "EOMES", "IL7R", "KLRB1", "ID2", "IL2RA", "THY1", "TBX21", "GATA3", "IL23R"), min.cutoff = "q9")
dev.off()

## Assign identity to clusters
new.cluster.ids <- c("T cell", "T cell", "T cell", "NK cell or ILC", "T cell", "NK cell or ILC", "T cell", "Indeterminate", "T or NK cell", "Indeterminate")
names(new.cluster.ids) <- levels(seurat.12.inter)
seurat.12.inter <- RenameIdents(seurat.12.inter, new.cluster.ids)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12inter-labeledUMAP.png")
DimPlot(seurat.12.inter, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
dev.off()

## Subset to cells with possible ILC identity (NK or ILC, NK or T)
seurat.12.ilc <- subset(seurat.12.inter, idents = c("NK cell or ILC", "T or NK cell"))

## Cluster the cells
DefaultAssay(seurat.12.ilc) <- "integrated"
seurat.12.ilc <- FindNeighbors(seurat.12.ilc, dims = 1:10)
seurat.12.ilc <- FindClusters(seurat.12.ilc, resolution = 0.5)

### Look at cluster IDs of the first 5 cells
head(Idents(seurat.12.ilc), 5)

## Run non-linear dimensional reduction (UMAP/tSNE)
seurat.12.ilc <- RunUMAP(seurat.12.ilc, dims = 1:10)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12ilc-initialUMAP.png")
DimPlot(seurat.12.ilc, reduction = "umap", label = TRUE)
dev.off()
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12ilc-initialUMAPcond.png")
DimPlot(seurat.12.ilc, reduction = "umap", group.by = "condition")
dev.off()

## Find markers that differentiate clusters from one another (to help assign identities)
DefaultAssay(seurat.12.ilc) <- "RNA"
seurat.12.ilc.markers <- FindAllMarkers(seurat.12.ilc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
manual.curate.12.ilc <- seurat.12.ilc.markers %>%
  group_by(cluster) %>%
  slice_max(n = 4, order_by = avg_log2FC)
write.table(manual.curate.12.ilc, file="manual.curate.12.ilc.txt", row.names=TRUE, col.names=TRUE)

## Visualize specific markers to identify cell types
### ILC1s: CD127/IL7R+, CD161/KLRB1+, ID2+, T-bet/TBX21+, CD294/PTGDR2-, CD117/KIT-
## ILC2s: CD127/IL7R+, CD161/KLRB1+, ID2+, CD294/PTGDR2+, GATA3+, CD117/KIT+-
### ILC3s: CD127/IL7R+, CD161/KLRB1+, ID2+, CD294/PTGDR-, CD117/KIT+, IL23R+
### NKs: CD56/NCAM1+, EOMES+
### Ts: CD3/CD3D+ (definite), TBX21+
### Monocytes: CD14+
### Neutrophils, Basophils, Eosinophils, DCs: CD11b+
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12ilc-initialmarkers.png", width= 700, height=480)
FeaturePlot(seurat.12.ilc, features = c("CD3D", "CD8A", "ITGAM", "CD14", "CD19", "NCAM1", "EOMES", "IL7R", "KLRB1", "ID2", "IL2RA", "THY1", "TBX21", "GATA3", "IL23R"), min.cutoff = "q9")
dev.off()

## Assign identity to clusters
## Cells do not have discernable levels of specific identifying markers; UMAP/ clustering did not separate T cells from NKs and ILCs by marker expression
## Manual marker ID
# new.cluster.ids <- c("")
# names(new.cluster.ids) <- levels(seurat.12.ilc)
# seurat.12.ilc <- RenameIdents(seurat.12.ilc, new.cluster.ids)
# png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12ilc-labeledUMAP.png")
# DimPlot(seurat.12.ilc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
# dev.off()

## Identify conserved markers
seurat.12.interc.markers <- FindConservedMarkers(seurat.12.inter, ident.1 = c("NK cell or ILC", "T or NK cell"), grouping.var = "condition", verbose = FALSE)
write.table(seurat.12.interc.markers, file="conserved.12.inter.txt", row.names=TRUE, col.names=TRUE)

## Visualize specific immune markers by condition and cell type
markers.to.plot <- c("CD3D", "CD8A", "ITGAM", "CD14", "CD19", "NCAM1", "EOMES", "IL7R", "KLRB1", "ID2", "IL2RA", "THY1", "TBX21", "GATA3", "IL23R")
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat*12inter-immdotplot.png")
DotPlot(seurat.12.inter, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8, split.by = "condition") +
  RotatedAxis()
dev.off()

## Visualize specific microbiome interaction genes by condition and cell type
markers.to.plot <- c("NOD2", "CARD9", "ATG16L1*", "IRGM")
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12inter-microdotplot.png")
DotPlot(seurat.12.inter, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8, split.by = "condition") +
  RotatedAxis()
dev.off()

seurat.12.inter$celltype <- Idents(seurat.12.inter)
plots <- VlnPlot(seurat.12.inter, features = c("NOD2", "CARD9", "ATG16L1", "IRGM"), split.by = "condition", group.by = "celltype",
    pt.size = 0, combine = FALSE)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12inter-vinmicrogenes.png")
wrap_plots(plots = plots, ncol = 1)
dev.off()

## Visualize specific microbiome interaction genes in cowplot format by condition and cell type
## All cells
theme_set(theme_cowplot())
Idents(seurat.12) <- "condition"
avg.seurat.12 <- as.data.frame(log1p(AverageExpression(seurat.12, verbose = FALSE)$RNA))
avg.seurat.12$gene <- rownames(avg.seurat.12)

genes.to.label = c("NOD2", "CARD9", "ATG16L1", "IRGM") 
p9 <- ggplot(avg.seurat.12, aes(healthy, UC)) + geom_point() + ggtitle("All CD45+ cells")
p9 <- LabelPoints(plot = p9, points = genes.to.label, repel = TRUE)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12-diffmicrogenes.png")
p9
dev.off()

## Subsetted cells
Idents(seurat.12.ilc) <- "condition"
avg.ilcs <- as.data.frame(log1p(AverageExpression(seurat.12.ilc, verbose = FALSE)$RNA))
avg.ilcs$gene <- rownames(avg.ilcs)

genes.to.label = c("NOD2", "CARD9", "ATG16L1", "IRGM") 
p9 <- ggplot(avg.ilcs, aes(healthy, UC)) + geom_point() + ggtitle("Selected immune cells")
p9 <- LabelPoints(plot = p9, points = genes.to.label, repel = TRUE)
png(file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/graphics/seurat12ilc-diffmicrogenes.png")
p9
dev.off()



## Save objects for future recall
saveRDS(seurat.12, file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat12-clustered.rds")
saveRDS(seurat.12.inter, file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat12inter-clustered.rds")
saveRDS(seurat.12.ilc, file = "C:/Users/Me/OneDrive - University of Nebraska at Omaha/Administrative/Documents/Senior Project/senior-capstone/rds/seurat12ilc-clustered.rds")

```

---
Metabolic pathway analysis
---

This section contains code for the KEGG metabolic pathway analysis sub-section. 

```{r Extract differential expression matrices}

# GSE 150050
## All-cluster matrix
im.15.markers <- FindAllMarkers(seurat.15, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
manual.curate.15 <- im.15.markers %>%
  group_by(cluster) %>%
  slice_max(n = 10, order_by = avg_log2FC)
im.15.markers <- im.15.markers %>% dplyr::select(p_val, avg_log2FC)
im.15.markers <- rownames_to_column(im.15.markers, "gene_id")
# Replace 0s
im.15.markers$p_val[im.15.markers$p_val=="0"]<-1.0e-302
im.15.markers <- im.15.markers[c("gene_id", "avg_log2FC", "p_val")]

### ILC versus other clusters matrix
im.ilc.markers <- FindMarkers(seurat.15.ilc, ident.1 = c("ILC1", "ILC2", "ILC3"), min.pct = 0.25)
im.ilc.markers <- im.ilc.markers %>% dplyr::select(p_val, avg_log2FC)
im.ilc.markers <- rownames_to_column(im.ilc.markers, "gene_id")
# Replace 0s
im.ilc.markers$p_val[im.ilc.markers$p_val=="0"]<-1.0e-302
im.ilc.markers <- im.ilc.markers[c("gene_id", "avg_log2FC", "p_val")]

### ILC subsets versus each other matrices
im.ilc1.markers <- FindMarkers(seurat.15.ilc, ident.1 = "ILC1", indent.2 = c("ILC2", "ILC3"), min.pct = 0.25)
im.ilc1.markers <- im.ilc1.markers %>% dplyr::select(p_val, avg_log2FC)
im.ilc1.markers <- rownames_to_column(im.ilc1.markers, "gene_id")
# Replace 0s
im.ilc1.markers$p_val[im.ilc1.markers$p_val=="0"]<-1.0e-302

im.ilc2.markers <- FindMarkers(seurat.15.ilc, ident.1 = "ILC2", indent.2 = c("ILC1", "ILC3"), min.pct = 0.25)
im.ilc2.markers <- im.ilc2.markers %>% dplyr::select(p_val, avg_log2FC)
im.ilc2.markers <- rownames_to_column(im.ilc2.markers, "gene_id")
# Replace 0s
im.ilc2.markers$p_val[im.ilc2.markers$p_val=="0"]<-1.0e-302

im.ilc3.markers <- FindMarkers(seurat.15.ilc, ident.1 = "ILC3", indent.2 = c("ILC1", "ILC2"), min.pct = 0.25)
im.ilc3.markers <- im.ilc3.markers %>% dplyr::select(p_val, avg_log2FC)
im.ilc3.markers <- rownames_to_column(im.ilc3.markers, "gene_id")
# Replace 0s
im.ilc3.markers$p_val[im.ilc3.markers$p_val=="0"]<-1.0e-302

# GSE 125527 
## Intestinal immune cells
## All-cluster matrix
int.12.markers <- FindAllMarkers(seurat.12.int, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
manual.curate.12.int <- int.12.markers %>%
  group_by(cluster) %>%
  slice_max(n = 10, order_by = avg_log2FC)
int.12.markers <- int.12.markers %>% dplyr::select(p_val, avg_log2FC)
int.12.markers <- rownames_to_column(int.12.markers, "gene_id")
# Replace 0s
int.12.markers$p_val[int.12.markers$p_val=="0"]<-1.0e-302

### ILC versus other clusters matrix
int.ilc.markers <- FindMarkers(seurat.12.int.ilc, ident.1 = "Contains ILCs", min.pct = 0.25)
int.ilc.markers <- int.ilc.markers %>% dplyr::select(p_val, avg_log2FC)
int.ilc.markers <- rownames_to_column(int.ilc.markers, "gene_id")
# Replace 0s
int.ilc.markers$p_val[int.ilc.markers$p_val=="0"]<-1.0e-302


# PBMCs
## All-cluster matrix
pbmc.12.markers <- FindAllMarkers(seurat.12.inter, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
manual.curate.12 <- pbmc.12.markers %>%
  group_by(cluster) %>%
  slice_max(n = 10, order_by = avg_log2FC)
pbmc.12.markers <- pbmc.12.markers %>% dplyr::select(p_val, avg_log2FC)
pbmc.12.markers <- rownames_to_column(pbmc.12.markers, "gene_id")
# Replace 0s
pbmc.12.markers$p_val[pbmc.12.markers$p_val=="0"]<-1.0e-302

### ILC versus other clusters matrix
inter.ilc.markers <- FindMarkers(seurat.12.inter, ident.1 = c("NK cell or ILC", "T or NK cell"), min.pct = 0.25)
inter.ilc.markers <- inter.ilc.markers %>% dplyr::select(p_val, avg_log2FC)
inter.ilc.markers <- rownames_to_column(inter.ilc.markers, "gene_id")
# Replace 0s
inter.ilc.markers$p_val[inter.ilc.markers$p_val=="0"]<-1.0e-302

## Differences by condition
### All cluster matrix, by condition
hu.12.markers <- FindMarkers(seurat.12, ident.1 = "healthy", ident.2 = "UC", verbose = FALSE)
hu.12.markers <- hu.12.markers %>% dplyr::select(p_val, avg_log2FC)
hu.12.markers <- rownames_to_column(hu.12.markers, "gene_id")
# Replace 0s
hu.12.markers$p_val[hu.12.markers$p_val=="0"]<-1.0e-302

### ILC matrix by condition
Idents(seurat.12.ilc) <- "condition"
hu.ilc.markers <- FindMarkers(seurat.12.ilc, ident.1 = "healthy", ident.2 = "UC",  min.pct = 0.25)
hu.ilc.markers <- hu.ilc.markers %>% dplyr::select(p_val, avg_log2FC)
hu.ilc.markers <- rownames_to_column(hu.ilc.markers, "gene_id")
# Replace 0s
hu.ilc.markers$p_val[hu.ilc.markers$p_val=="0"]<-1.0e-302


```

```{r Querying against KEGG}

## KEGG
## GSE 150050
im.15.path <- run_pathfindR(im.15.markers, output_dir = "./pathfindR/seurat15_all")
im.ilc.path <- run_pathfindR(im.ilc.markers, output_dir = "./pathfindR/seurat15_ilc")
im.ilc1.path <- run_pathfindR(im.ilc1.markers, output_dir = "./pathfindR/seurat15_ilc1")
im.ilc2.path <- run_pathfindR(im.ilc2.markers, output_dir = "./pathfindR/seurat15_ilc2")
im.ilc3.path <- run_pathfindR(im.ilc3.markers, output_dir = "./pathfindR/seurat15_ilc2")

## GSE 125527
## Intestinal immune cells
int.12.path <- run_pathfindR(int.12.markers, output_dir = "./pathfindR/seurat12int_all")
int.ilc.path <- run_pathfindR(int.ilc.markers, output_dir = "./pathfindR/seurat12int_ilc")

## PBMCs 
pbmc.12.path <- run_pathfindR(pbmc.12.markers, output_dir = "./pathfindR/seurat12_all")
inter.ilc.path <- run_pathfindR(inter.ilc.markers, output_dir = "./pathfindR/seurat12_ilc")

## PBMCs by condition
hu.12.path <- run_pathfindR(hu.12.markers, output_dir = "./pathfindR/seurat12_conditions")
hu.ilc.path <- run_pathfindR(hu.ilc.markers, output_dir = "./pathfindR/seurat12_ilc_conditions")

## GO-BP
## GSE 150050
im.15.path <- run_pathfindR(im.15.markers, output_dir = "./pathfindR/seurat15_all", gene_sets = "GO-BP")
im.ilc.path <- run_pathfindR(im.ilc.markers, output_dir = "./pathfindR/seurat15_ilc", gene_sets = "GO-BP")
im.ilc1.path <- run_pathfindR(im.ilc1.markers, output_dir = "./pathfindR/seurat15_ilc1", gene_sets = "GO-BP")
im.ilc2.path <- run_pathfindR(im.ilc2.markers, output_dir = "./pathfindR/seurat15_ilc2", gene_sets = "GO-BP")
im.ilc3.path <- run_pathfindR(im.ilc3.markers, output_dir = "./pathfindR/seurat15_ilc2", gene_sets = "GO-BP")

## GSE 125527
## Intestinal immune cells
int.12.path <- run_pathfindR(int.12.markers, output_dir = "./pathfindR/seurat12int_all", gene_sets = "GO-BP")
int.ilc.path <- run_pathfindR(int.ilc.markers, output_dir = "./pathfindR/seurat12int_ilc", gene_sets = "GO-BP")

## PBMCs 
pbmc.12.path <- run_pathfindR(pbmc.12.markers, output_dir = "./pathfindR/seurat12_all", gene_sets = "GO-BP")
inter.ilc.path <- run_pathfindR(inter.ilc.markers, output_dir = "./pathfindR/seurat12_ilc", gene_sets = "GO-BP")

## PBMCs by condition
hu.12.path <- run_pathfindR(hu.12.markers, output_dir = "./pathfindR/seurat12_conditions", gene_sets = "GO-BP")
hu.ilc.path <- run_pathfindR(hu.ilc.markers, output_dir = "./pathfindR/seurat12_ilc_conditions", gene_sets = "GO-BP")


# Visualize metabolic genes using dot plot
markers.to.plot <- c("NOD2", "CARD9", "ATG16L1", "IRGM")
png(file = "~/graphics/~.png")
DotPlot(seurat, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8) +
  RotatedAxis()

```