processing.Rmd

---
title: "Processing and Demultiplexing"
output: html_notebook
---

```{r}
library(Seurat)
library(cowplot)
library(dplyr)
library(tidyr)
library(viridis)
library(KernSmooth)
library(ggplot2)
library(reshape2)
library(Rtsne)
library(ShortRead)
library(stringdist)
library(deMULTIplex)
```

# Load in data
Because of multiplexing, the samples are distributed a little weird.

Here are the directories containing each sample:
P3_1: P3_P7_1 and P3_P7_2 - (1)
P3_2: P3_2 - (2)
P7_1: P7_1 - (3)
P7_2: P3_P7_1 and P3_P7_2 - (1)
P14_1: P14_1 - (4)
P14_2: P14_2 - (2)

Numbers in parantheses represent processing day

```{r}
P3_P7_1_mat <- Read10X(data.dir = "../data/cellranger/P3_P7_1/filtered_feature_bc_matrix/")
P3_P7_2_mat <- Read10X(data.dir = "../data/cellranger/P3_P7_2/filtered_feature_bc_matrix/")

P3_2_mat <- Read10X(data.dir = "../data/cellranger/P3_2/filtered_feature_bc_matrix/")

P7_1_mat <- Read10X(data.dir = "../data/cellranger/P7_1/filtered_feature_bc_matrix/")

P14_1_mat <- Read10X(data.dir = "../data/cellranger/P14_1/filtered_feature_bc_matrix/")
P14_2_mat <- Read10X(data.dir = "../data/cellranger/P14_2/filtered_feature_bc_matrix/")
```

# Quick seurat processing
## P3_P7_1
```{r}
seurat <- CreateSeuratObject(P3_P7_1_mat, min.cells = 3, 
                             min.features = 200,
                             project = "P3_P7_1")
```

```{r}
mito.features <- grep(pattern = "^mt-", x = rownames(seurat), value = TRUE)
percent.mito <- Matrix::colSums(GetAssayData(seurat, slot = 'counts')[mito.features, ]) / Matrix::colSums(GetAssayData(seurat, slot = 'counts'))

# The [[ operator can add columns to object metadata, and is a great place to stash QC stats
seurat[['percent.mito']] <- percent.mito
VlnPlot(seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"), ncol = 3)
```

```{r}
hist(seurat$percent.mito, breaks=50)
```

```{r}
seurat <- subset(seurat, 
                 subset = nFeature_RNA > 200 & 
                   percent.mito < 0.2)
```

```{r}
P3_P7_1_seurat <- seurat
P3_P7_1_keep <- colnames(P3_P7_1_seurat)
```

## P3_P7_2
```{r}
seurat <- CreateSeuratObject(P3_P7_2_mat, min.cells = 3, 
                             min.features = 200,
                             project = "P3_P7_2")
```

```{r}
mito.features <- grep(pattern = "^mt-", x = rownames(seurat), value = TRUE)
percent.mito <- Matrix::colSums(GetAssayData(seurat, slot = 'counts')[mito.features, ]) / Matrix::colSums(GetAssayData(seurat, slot = 'counts'))

# The [[ operator can add columns to object metadata, and is a great place to stash QC stats
seurat[['percent.mito']] <- percent.mito
VlnPlot(seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"), ncol = 3)
```

```{r}
hist(seurat$percent.mito, breaks=50)
```

```{r}
seurat <- subset(seurat, 
                 subset = nFeature_RNA > 200 & 
                   percent.mito < 0.2)
```

```{r}
P3_P7_2_seurat <- seurat
P3_P7_2_keep <- colnames(P3_P7_2_seurat)
```

## P3_2
```{r}
seurat <- CreateSeuratObject(P3_2_mat, min.cells = 3, 
                             min.features = 200,
                             project = "P3_2")
```

```{r}
mito.features <- grep(pattern = "^mt-", x = rownames(seurat), value = TRUE)
percent.mito <- Matrix::colSums(GetAssayData(seurat, slot = 'counts')[mito.features, ]) / Matrix::colSums(GetAssayData(seurat, slot = 'counts'))

# The [[ operator can add columns to object metadata, and is a great place to stash QC stats
seurat[['percent.mito']] <- percent.mito
VlnPlot(seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"), ncol = 3)
```

```{r}
hist(seurat$percent.mito, breaks=50)
```

```{r}
seurat <- subset(seurat, 
                 subset = nFeature_RNA > 200 & 
                   percent.mito < 0.2)
```

```{r}
P3_2_seurat <- seurat
P3_2_keep <- colnames(P3_2_seurat)
```

## P7_1
```{r}
seurat <- CreateSeuratObject(P7_1_mat, min.cells = 3, 
                             min.features = 200,
                             project = "P7_1")
```

```{r}
mito.features <- grep(pattern = "^mt-", x = rownames(seurat), value = TRUE)
percent.mito <- Matrix::colSums(GetAssayData(seurat, slot = 'counts')[mito.features, ]) / Matrix::colSums(GetAssayData(seurat, slot = 'counts'))

# The [[ operator can add columns to object metadata, and is a great place to stash QC stats
seurat[['percent.mito']] <- percent.mito
VlnPlot(seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"), ncol = 3)
```

```{r}
hist(seurat$percent.mito, breaks=50)
```

```{r}
seurat <- subset(seurat, 
                 subset = nFeature_RNA > 200 & 
                   percent.mito < 0.2)
```

```{r}
P7_1_seurat <- seurat
P7_1_keep <- colnames(P7_1_seurat)
```

## P14_1
```{r}
seurat <- CreateSeuratObject(P14_1_mat, min.cells = 3, 
                             min.features = 200,
                             project = "P14_1")
```

```{r}
mito.features <- grep(pattern = "^mt-", x = rownames(seurat), value = TRUE)
percent.mito <- Matrix::colSums(GetAssayData(seurat, slot = 'counts')[mito.features, ]) / Matrix::colSums(GetAssayData(seurat, slot = 'counts'))

# The [[ operator can add columns to object metadata, and is a great place to stash QC stats
seurat[['percent.mito']] <- percent.mito
VlnPlot(seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"), ncol = 3)
```

```{r}
hist(seurat$percent.mito, breaks=50)
```


```{r}
seurat <- subset(seurat, 
                 subset = nFeature_RNA > 200 & 
                   percent.mito < 0.25)
```

```{r}
P14_1_seurat <- seurat
P14_1_keep <- colnames(P14_1_seurat)
```

## P14_2
```{r}
seurat <- CreateSeuratObject(P14_2_mat, min.cells = 3, 
                             min.features = 200,
                             project = "P14_2")
```

```{r}
mito.features <- grep(pattern = "^mt-", x = rownames(seurat), value = TRUE)
percent.mito <- Matrix::colSums(GetAssayData(seurat, slot = 'counts')[mito.features, ]) / Matrix::colSums(GetAssayData(seurat, slot = 'counts'))

# The [[ operator can add columns to object metadata, and is a great place to stash QC stats
seurat[['percent.mito']] <- percent.mito
VlnPlot(seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mito"), ncol = 3)
```

```{r}
hist(seurat$percent.mito, breaks=50)
```

```{r}
seurat <- subset(seurat, 
                 subset = nFeature_RNA > 200 & 
                   percent.mito < 0.25)
```

```{r}
P14_2_seurat <- seurat
P14_2_keep <- colnames(P14_2_seurat)
```

# MULTI-seq annotations
```{r}
bar.ref <- read.csv("../data/LMOlist.csv", stringsAsFactors = F,
                    header=F)$V1
```

## Quantify barcodes from fastq files
```{r}
#P3_P7_1
readTable <- MULTIseq.preProcess(R1 = "../fastq/P3_P7_1_Barcode/LungP3P7_1_Barcode_R1_paired.fastq.gz",
                                 R2 = "../fastq/P3_P7_1_Barcode/LungP3P7_1_Barcode_R2_paired.fastq.gz",
                                 cellIDs = P3_P7_1_keep,
                                 cell=c(1,16),
                                 umi=c(17,28),
                                 tag=c(1,8))
P3_P7_1_bar_table <- MULTIseq.align(readTable, P3_P7_1_keep, bar.ref[13:24])
write.csv(P3_P7_1_bar_table, file="../data/P3_P7_1_barcode_counts.csv", quote=F)

#P3_P7_2
readTable <- MULTIseq.preProcess(R1 = "../fastq/P3_P7_2_Barcode/LungP3P7_2_Barcode_R1_paired.fastq.gz",
                                 R2 = "../fastq/P3_P7_2_Barcode/LungP3P7_2_Barcode_R2_paired.fastq.gz",
                                 cellIDs = P3_P7_2_keep,
                                 cell=c(1,16),
                                 umi=c(17,28),
                                 tag=c(1,8))
P3_P7_2_bar_table <- MULTIseq.align(readTable, P3_P7_2_keep, bar.ref[13:24])
write.csv(P3_P7_2_bar_table, file="../data/P3_P7_2_barcode_counts.csv", quote=F)

#P3_2
readTable <- MULTIseq.preProcess(R1 = "../fastq/P3_2_Barcode/Lung_P3_2_Barcode_R1_paired.fastq.gz",
                                 R2 = "../fastq/P3_2_Barcode/Lung_P3_2_Barcode_R2_paired.fastq.gz",
                                 cellIDs = P3_2_keep,
                                 cell=c(1,16),
                                 umi=c(17,28),
                                 tag=c(1,8))
P3_2_bar_table <- MULTIseq.align(readTable, P3_2_keep, bar.ref[c(1,2,3,4,5,8)])
write.csv(P3_2_bar_table, file="../data/P3_2_barcode_counts.csv", quote=F)

#P7_1
readTable <- MULTIseq.preProcess(R1 = "../fastq/P7_1/Lung_P7_Barcode_R1_paired.fastq.gz",
                                 R2 = "../fastq/P7_1/Lung_P7_Barcode_R2_paired.fastq.gz",
                                 cellIDs = P7_1_keep,
                                 cell=c(1,16),
                                 umi=c(17,28),
                                 tag=c(1,8))
P7_1_bar_table <- MULTIseq.align(readTable, P7_1_keep, bar.ref[1:6])
write.csv(P7_1_bar_table, file="../data/P7_1_barcode_counts.csv", quote=F)

#P14_1
readTable <- MULTIseq.preProcess(R1 = "../fastq/P14_1_Barcode/Lung_P14_Barcode_R1_paired.fastq.gz",
                                 R2 = "../fastq/P14_1_Barcode/Lung_P14_Barcode_R2_paired.fastq.gz",
                                 cellIDs = P14_1_keep,
                                 cell=c(1,16),
                                 umi=c(17,28),
                                 tag=c(1,8))
P14_1_bar_table <- MULTIseq.align(readTable, P14_1_keep, bar.ref[25:30])
write.csv(P14_1_bar_table, file="../data/P14_1_barcode_counts.csv", quote=F)

#P14_2
readTable <- MULTIseq.preProcess(R1 = "../fastq/P14_2_Barcode/Lung_P14_2_Barcode_R1_paired.fastq.gz",
                                 R2 = "../fastq/P14_2_Barcode/Lung_P14_2_Barcode_R2_paired.fastq.gz",
                                 cellIDs = P14_2_keep,
                                 cell=c(1,16),
                                 umi=c(17,28),
                                 tag=c(1,8))
P14_2_bar_table <- MULTIseq.align(readTable, P14_2_keep, bar.ref[7:12])
write.csv(P14_2_bar_table, file="../data/P14_2_barcode_counts.csv", quote=F)
```

# Annotate samples
## P3_P7_1
```{r}
tsne <- barTSNE(P3_P7_1_bar_table[,1:12])
```

```{r}
plot(tsne$TSNE1, tsne$TSNE2, pch=20)
```

```{r}
temp <- P3_P7_1_bar_table
temp$TSNE1 <- tsne$TSNE1
temp$TSNE2 <- tsne$TSNE2
write.csv(temp, file="../output/P3_P7_1_barcode_tsne.csv", quote=F)
```

```{r}
temp <- as.matrix(tsne[3:14])
temp[temp<0] <- 0
temp <- as.data.frame(temp)
temp$TSNE1 <- tsne$TSNE1
temp$TSNE2 <- tsne$TSNE2


for (LMO in colnames(tsne)[3:14]){
  png(filename = paste0("../figs/LMOs/P3_P7_1/",LMO,".png"), width = 700, height = 600)
  p <- ggplot(data=temp, aes_string(x="TSNE1",y="TSNE2",color=LMO)) + geom_point() +
    scale_color_gradient(low="lightgrey",high="red") +
    theme_void()
  print(p)
  dev.off()
}
```

```{r}
bar_filtered <- P3_P7_1_bar_table[,1:12]
#bar_filtered <- bar_filtered[,-bad_barcodes]
```

```{r}
bar.table_sweep.list <- list()
n <- 0
for (q in seq(0.01, 0.99, by=0.02)) {
  print(q)
  n <- n + 1
  bar.table_sweep.list[[n]] <- classifyCells(bar_filtered, q=q)
  names(bar.table_sweep.list)[n] <- paste("q=",q,sep="")
}
```

```{r}
threshold.results1 <- findThresh(call.list=bar.table_sweep.list)
ggplot(data=threshold.results1$res, aes(x=q, y=Proportion, color=Subset)) + geom_line() + theme(legend.position = "none") +
  geom_vline(xintercept=threshold.results1$extrema, lty=2) + scale_color_manual(values=c("red","black","blue"))
```

```{r}
round1.calls <- classifyCells(bar_filtered, 
                              q=findQ(threshold.results1$res, 
                                      threshold.results1$extrema))
```

View classifications
```{r}
tsne$Classification <- "Singlet"
tsne$Classification[which(round1.calls[rownames(tsne)]=="Doublet")] <- "Doublet"
tsne$Classification[which(round1.calls[rownames(tsne)]=="Negative")] <- "Negative"
```

```{r}
table(tsne$Classification)
```


```{r}
tsne_classification <- ggplot(tsne, aes(x=TSNE1, y=TSNE2)) +
  geom_point(size=0.25, aes(color=Classification)) +
  theme_void()
ggsave(tsne_classification, file="../figs/P3_P7_1_classification_tsne.png",
       dpi=600, width=8, height=5.5)
tsne_classification
```

Reclassify negative cells
```{r}
bar_full <- P3_P7_1_bar_table[,1:12]
reclass.cells <- findReclassCells(bar_full,
                                  unique(names(round1.calls)[which(round1.calls=="Negative")]))
reclass.res <- rescueCells(bar_full, round1.calls[rownames(bar_full)], reclass.cells)
```

```{r}
ggplot(reclass.res[-1, ], aes(x=ClassStability, y=MatchRate_mean)) + 
    geom_point() + xlim(c(nrow(reclass.res)-1,1)) + 
    ylim(c(0,1.05)) +
    geom_errorbar(aes(ymin=MatchRate_mean-MatchRate_sd, ymax=MatchRate_mean+MatchRate_sd), width=.1) +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1], color="red") +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1]+3*reclass.res$MatchRate_sd[1], color="red",lty=2) +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1]-3*reclass.res$MatchRate_sd[1], color="red",lty=2)
```

```{r}
final.calls <- round1.calls
rescue.ind <- which(reclass.cells$ClassStability >= 8) ## Note: Value will be dataset-specific
final.calls[rownames(reclass.cells)[rescue.ind]] <- reclass.cells$Reclassification[rescue.ind]
```

Recheck tSNE
```{r}
tsne$Classification <- "Singlet"
tsne$Classification[which(final.calls[rownames(tsne)]=="Doublet")] <- "Doublet"
tsne$Classification[which(final.calls[rownames(tsne)]=="Negative")] <- "Negative"
```

```{r}
table(tsne$Classification)
```

```{r}
tsne_classification <- ggplot(tsne, aes(x=TSNE1, y=TSNE2)) +
  geom_point(size=0.25, aes(color=Classification)) +
  theme_void()
ggsave(tsne_classification, file="../figs/P3_P7_1_classification_tsne.png",
       dpi=600, width=8, height=5.5)
tsne_classification
```

Plot annotations on tsne plot
```{r}
samples <- unique(final.calls)
plotSampleTSNE <- function(sample){
  data <- tsne
  data$Sample <- "Other"
  data$Sample[which(final.calls[rownames(tsne)]==sample)] <- sample
  sample_plot <- ggplot(data, aes(x=TSNE1, y=TSNE2)) +
    geom_point(size=0.25, alpha=0.5, aes(color=Sample)) +
    scale_color_manual(values=c("red", "lightgrey")) +
    theme_void()
  ggsave(sample_plot, file=paste0("../figs/LMOs/P3_P7_1/Classifications/", 
         sample,".png"), width=5, height=3.2)
}
```

```{r}
lapply(samples, plotSampleTSNE)
```

```{r}
P3_P7_1_calls <- final.calls
```


## P3_P7_2

```{r}
tsne <- barTSNE(P3_P7_2_bar_table[,1:12])
```

```{r}
plot(tsne$TSNE1, tsne$TSNE2, pch=20)
```

```{r}
temp <- P3_P7_2_bar_table
temp$TSNE1 <- tsne$TSNE1
temp$TSNE2 <- tsne$TSNE2
write.csv(temp, file="../output/P3_P7_2_barcode_tsne.csv", quote=F)
```

```{r}
temp <- as.matrix(tsne[3:14])
temp[temp<0] <- 0
temp <- as.data.frame(temp)
temp$TSNE1 <- tsne$TSNE1
temp$TSNE2 <- tsne$TSNE2


for (LMO in colnames(tsne)[3:14]){
  png(filename = paste0("../figs/LMOs/P3_P7_2/",LMO,".png"), width = 700, height = 600)
  p <- ggplot(data=temp, aes_string(x="TSNE1",y="TSNE2",color=LMO)) + geom_point() +
    scale_color_gradient(low="lightgrey",high="red") +
    theme_void()
  print(p)
  dev.off()
}
```

```{r}
bar_filtered <- P3_P7_2_bar_table[,1:12]
#bar_filtered <- bar_filtered[,-bad_barcodes]
```

```{r}
bar.table_sweep.list <- list()
n <- 0
for (q in seq(0.01, 0.99, by=0.02)) {
  print(q)
  n <- n + 1
  bar.table_sweep.list[[n]] <- classifyCells(bar_filtered, q=q)
  names(bar.table_sweep.list)[n] <- paste("q=",q,sep="")
}
```

```{r}
threshold.results1 <- findThresh(call.list=bar.table_sweep.list)
ggplot(data=threshold.results1$res, aes(x=q, y=Proportion, color=Subset)) + geom_line() + theme(legend.position = "none") +
  geom_vline(xintercept=threshold.results1$extrema, lty=2) + scale_color_manual(values=c("red","black","blue"))
```

```{r}
round1.calls <- classifyCells(bar_filtered, 
                              q=findQ(threshold.results1$res, 
                                      threshold.results1$extrema))
```

View classifications
```{r}
tsne$Classification <- "Singlet"
tsne$Classification[which(round1.calls[rownames(tsne)]=="Doublet")] <- "Doublet"
tsne$Classification[which(round1.calls[rownames(tsne)]=="Negative")] <- "Negative"
```

```{r}
table(tsne$Classification)
```


```{r}
tsne_classification <- ggplot(tsne, aes(x=TSNE1, y=TSNE2)) +
  geom_point(size=0.25, aes(color=Classification)) +
  theme_void()
ggsave(tsne_classification, file="../figs/P3_P7_2_classification_tsne.png",
       dpi=600, width=8, height=5.5)
tsne_classification
```

Reclassify negative cells
```{r}
bar_full <- P3_P7_2_bar_table[,1:12]
reclass.cells <- findReclassCells(bar_full,
                                  unique(names(round1.calls)[which(round1.calls=="Negative")]))
reclass.res <- rescueCells(bar_full, round1.calls[rownames(bar_full)], reclass.cells)
```

```{r}
ggplot(reclass.res[-1, ], aes(x=ClassStability, y=MatchRate_mean)) + 
    geom_point() + xlim(c(nrow(reclass.res)-1,1)) + 
    ylim(c(0,1.05)) +
    geom_errorbar(aes(ymin=MatchRate_mean-MatchRate_sd, ymax=MatchRate_mean+MatchRate_sd), width=.1) +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1], color="red") +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1]+3*reclass.res$MatchRate_sd[1], color="red",lty=2) +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1]-3*reclass.res$MatchRate_sd[1], color="red",lty=2)
```

```{r}
final.calls <- round1.calls
#Class stability seems low, we won't re-classify
#rescue.ind <- which(reclass.cells$ClassStability >= 8) ## Note: Value will be dataset-specific
#final.calls[rownames(reclass.cells)[rescue.ind]] <- reclass.cells$Reclassification[rescue.ind]
```

Recheck tSNE
```{r}
tsne$Classification <- "Singlet"
tsne$Classification[which(final.calls[rownames(tsne)]=="Doublet")] <- "Doublet"
tsne$Classification[which(final.calls[rownames(tsne)]=="Negative")] <- "Negative"
```

```{r}
table(tsne$Classification)
```

```{r}
tsne_classification <- ggplot(tsne, aes(x=TSNE1, y=TSNE2)) +
  geom_point(size=0.25, aes(color=Classification)) +
  theme_void()
ggsave(tsne_classification, file="../figs/P3_P7_2_classification_tsne.png",
       dpi=600, width=8, height=5.5)
tsne_classification
```

Plot annotations on tsne plot
```{r}
samples <- unique(final.calls)
plotSampleTSNE <- function(sample){
  data <- tsne
  data$Sample <- "Other"
  data$Sample[which(final.calls[rownames(tsne)]==sample)] <- sample
  sample_plot <- ggplot(data, aes(x=TSNE1, y=TSNE2)) +
    geom_point(size=0.25, alpha=0.5, aes(color=Sample)) +
    scale_color_manual(values=c("red", "lightgrey")) +
    theme_void()
  ggsave(sample_plot, file=paste0("../figs/LMOs/P3_P7_2/Classifications/", 
         sample,".png"), width=5, height=3.2)
}
```

```{r}
lapply(samples, plotSampleTSNE)
```

```{r}
P3_P7_2_calls <- final.calls
```

## P3_2
```{r}
tsne <- barTSNE(P3_2_bar_table[,1:6])
```

```{r}
plot(tsne$TSNE1, tsne$TSNE2, pch=20)
```

```{r}
temp <- P3_2_bar_table
temp$TSNE1 <- tsne$TSNE1
temp$TSNE2 <- tsne$TSNE2
write.csv(temp, file="../output/P3_2_barcode_tsne.csv", quote=F)
```

```{r}
temp <- as.matrix(tsne[3:8])
temp[temp<0] <- 0
temp <- as.data.frame(temp)
temp$TSNE1 <- tsne$TSNE1
temp$TSNE2 <- tsne$TSNE2


for (LMO in colnames(tsne)[3:8]){
  png(filename = paste0("../figs/LMOs/P3_2/",LMO,".png"), width = 700, height = 600)
  p <- ggplot(data=temp, aes_string(x="TSNE1",y="TSNE2",color=LMO)) + geom_point() +
    scale_color_gradient(low="lightgrey",high="red") +
    theme_void()
  print(p)
  dev.off()
}
```

```{r}
bar_filtered <- P3_2_bar_table[,1:6]
#bar_filtered <- bar_filtered[,-bad_barcodes]
```

```{r}
bar.table_sweep.list <- list()
n <- 0
for (q in seq(0.01, 0.99, by=0.02)) {
  print(q)
  n <- n + 1
  bar.table_sweep.list[[n]] <- classifyCells(bar_filtered, q=q)
  names(bar.table_sweep.list)[n] <- paste("q=",q,sep="")
}
```

```{r}
threshold.results1 <- findThresh(call.list=bar.table_sweep.list)
ggplot(data=threshold.results1$res, aes(x=q, y=Proportion, color=Subset)) + geom_line() + theme(legend.position = "none") +
  geom_vline(xintercept=threshold.results1$extrema, lty=2) + scale_color_manual(values=c("red","black","blue"))
```

```{r}
round1.calls <- classifyCells(bar_filtered, 
                              q=findQ(threshold.results1$res, 
                                      threshold.results1$extrema))
```

View classifications
```{r}
tsne$Classification <- "Singlet"
tsne$Classification[which(round1.calls[rownames(tsne)]=="Doublet")] <- "Doublet"
tsne$Classification[which(round1.calls[rownames(tsne)]=="Negative")] <- "Negative"
```

```{r}
table(tsne$Classification)
```


```{r}
tsne_classification <- ggplot(tsne, aes(x=TSNE1, y=TSNE2)) +
  geom_point(size=0.25, aes(color=Classification)) +
  theme_void()
ggsave(tsne_classification, file="../figs/P3_2_classification_tsne.png",
       dpi=600, width=8, height=5.5)
tsne_classification
```

Reclassify negative cells
```{r}
bar_full <- P3_2_bar_table[,1:6]
reclass.cells <- findReclassCells(bar_full,
                                  unique(names(round1.calls)[which(round1.calls=="Negative")]))
reclass.res <- rescueCells(bar_full, round1.calls[rownames(bar_full)], reclass.cells)
```

```{r}
ggplot(reclass.res[-1, ], aes(x=ClassStability, y=MatchRate_mean)) + 
    geom_point() + xlim(c(nrow(reclass.res)-1,1)) + 
    ylim(c(0,1.05)) +
    geom_errorbar(aes(ymin=MatchRate_mean-MatchRate_sd, ymax=MatchRate_mean+MatchRate_sd), width=.1) +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1], color="red") +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1]+3*reclass.res$MatchRate_sd[1], color="red",lty=2) +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1]-3*reclass.res$MatchRate_sd[1], color="red",lty=2)
```

```{r}
final.calls <- round1.calls
rescue.ind <- which(reclass.cells$ClassStability >= 13) ## Note: Value will be dataset-specific
final.calls[rownames(reclass.cells)[rescue.ind]] <- reclass.cells$Reclassification[rescue.ind]
```

Recheck tSNE
```{r}
tsne$Classification <- "Singlet"
tsne$Classification[which(final.calls[rownames(tsne)]=="Doublet")] <- "Doublet"
tsne$Classification[which(final.calls[rownames(tsne)]=="Negative")] <- "Negative"
```

```{r}
table(tsne$Classification)
```

```{r}
tsne_classification <- ggplot(tsne, aes(x=TSNE1, y=TSNE2)) +
  geom_point(size=0.25, aes(color=Classification)) +
  theme_void()
ggsave(tsne_classification, file="../figs/P3_2_classification_tsne.png",
       dpi=600, width=8, height=5.5)
tsne_classification
```

Plot annotations on tsne plot
```{r}
samples <- unique(final.calls)
plotSampleTSNE <- function(sample){
  data <- tsne
  data$Sample <- "Other"
  data$Sample[which(final.calls[rownames(tsne)]==sample)] <- sample
  sample_plot <- ggplot(data, aes(x=TSNE1, y=TSNE2)) +
    geom_point(size=0.25, alpha=0.5, aes(color=Sample)) +
    scale_color_manual(values=c("red", "lightgrey")) +
    theme_void()
  ggsave(sample_plot, file=paste0("../figs/LMOs/P3_2/Classifications/", 
         sample,".png"), width=5, height=3.2)
}
```

```{r}
lapply(samples, plotSampleTSNE)
```

```{r}
P3_2_calls <- final.calls
```

## P7_1
```{r}
tsne <- barTSNE(P7_1_bar_table[,1:6])
```

```{r}
plot(tsne$TSNE1, tsne$TSNE2, pch=20)
```

```{r}
temp <- P7_1_bar_table
temp$TSNE1 <- tsne$TSNE1
temp$TSNE2 <- tsne$TSNE2
write.csv(temp, file="../output/P7_1_barcode_tsne.csv", quote=F)
```

```{r}
temp <- as.matrix(tsne[3:8])
temp[temp<0] <- 0
temp <- as.data.frame(temp)
temp$TSNE1 <- tsne$TSNE1
temp$TSNE2 <- tsne$TSNE2


for (LMO in colnames(tsne)[3:8]){
  png(filename = paste0("../figs/LMOs/P7_1/",LMO,".png"), width = 700, height = 600)
  p <- ggplot(data=temp, aes_string(x="TSNE1",y="TSNE2",color=LMO)) + geom_point() +
    scale_color_gradient(low="lightgrey",high="red") +
    theme_void()
  print(p)
  dev.off()
}
```

```{r}
bar_filtered <- P7_1_bar_table[,1:6]
#bar_filtered <- bar_filtered[,-bad_barcodes]
```

```{r}
bar.table_sweep.list <- list()
n <- 0
for (q in seq(0.01, 0.99, by=0.02)) {
  print(q)
  n <- n + 1
  bar.table_sweep.list[[n]] <- classifyCells(bar_filtered, q=q)
  names(bar.table_sweep.list)[n] <- paste("q=",q,sep="")
}
```

```{r}
threshold.results1 <- findThresh(call.list=bar.table_sweep.list)
ggplot(data=threshold.results1$res, aes(x=q, y=Proportion, color=Subset)) + geom_line() + theme(legend.position = "none") +
  geom_vline(xintercept=threshold.results1$extrema, lty=2) + scale_color_manual(values=c("red","black","blue"))
```

```{r}
round1.calls <- classifyCells(bar_filtered, 
                              q=findQ(threshold.results1$res, 
                                      threshold.results1$extrema))
```

View classifications
```{r}
tsne$Classification <- "Singlet"
tsne$Classification[which(round1.calls[rownames(tsne)]=="Doublet")] <- "Doublet"
tsne$Classification[which(round1.calls[rownames(tsne)]=="Negative")] <- "Negative"
```

```{r}
table(tsne$Classification)
```


```{r}
tsne_classification <- ggplot(tsne, aes(x=TSNE1, y=TSNE2)) +
  geom_point(size=0.25, aes(color=Classification)) +
  theme_void()
ggsave(tsne_classification, file="../figs/P7_1_classification_tsne.png",
       dpi=600, width=8, height=5.5)
tsne_classification
```

Reclassify negative cells
```{r}
bar_full <- P7_1_bar_table[,1:6]
reclass.cells <- findReclassCells(bar_full,
                                  unique(names(round1.calls)[which(round1.calls=="Negative")]))
reclass.res <- rescueCells(bar_full, round1.calls[rownames(bar_full)], reclass.cells)
```

```{r}
ggplot(reclass.res[-1, ], aes(x=ClassStability, y=MatchRate_mean)) + 
    geom_point() + xlim(c(nrow(reclass.res)-1,1)) + 
    ylim(c(0,1.05)) +
    geom_errorbar(aes(ymin=MatchRate_mean-MatchRate_sd, ymax=MatchRate_mean+MatchRate_sd), width=.1) +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1], color="red") +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1]+3*reclass.res$MatchRate_sd[1], color="red",lty=2) +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1]-3*reclass.res$MatchRate_sd[1], color="red",lty=2)
```

Re-classification doesn't seem promising
```{r}
final.calls <- round1.calls
```

Recheck tSNE
```{r}
tsne$Classification <- "Singlet"
tsne$Classification[which(final.calls[rownames(tsne)]=="Doublet")] <- "Doublet"
tsne$Classification[which(final.calls[rownames(tsne)]=="Negative")] <- "Negative"
```

```{r}
table(tsne$Classification)
```

```{r}
tsne_classification <- ggplot(tsne, aes(x=TSNE1, y=TSNE2)) +
  geom_point(size=0.25, aes(color=Classification)) +
  theme_void()
ggsave(tsne_classification, file="../figs/P7_1_classification_tsne.png",
       dpi=600, width=8, height=5.5)
tsne_classification
```

Plot annotations on tsne plot
```{r}
samples <- unique(final.calls)
plotSampleTSNE <- function(sample){
  data <- tsne
  data$Sample <- "Other"
  data$Sample[which(final.calls[rownames(tsne)]==sample)] <- sample
  sample_plot <- ggplot(data, aes(x=TSNE1, y=TSNE2)) +
    geom_point(size=0.25, alpha=0.5, aes(color=Sample)) +
    scale_color_manual(values=c("red", "lightgrey")) +
    theme_void()
  ggsave(sample_plot, file=paste0("../figs/LMOs/P7_1/Classifications/", 
         sample,".png"), width=5, height=3.2)
}
```

```{r}
lapply(samples, plotSampleTSNE)
```

```{r}
P7_1_calls <- final.calls
```

## P14_1
```{r}
tsne <- barTSNE(P14_1_bar_table[,1:6])
```

```{r}
plot(tsne$TSNE1, tsne$TSNE2, pch=20)
```

```{r}
temp <- P14_1_bar_table
temp$TSNE1 <- tsne$TSNE1
temp$TSNE2 <- tsne$TSNE2
write.csv(temp, file="../output/P14_1_barcode_tsne.csv", quote=F)
```

```{r}
temp <- as.matrix(tsne[3:8])
temp[temp<0] <- 0
temp <- as.data.frame(temp)
temp$TSNE1 <- tsne$TSNE1
temp$TSNE2 <- tsne$TSNE2


for (LMO in colnames(tsne)[3:8]){
  png(filename = paste0("../figs/LMOs/P14_1/",LMO,".png"), width = 700, height = 600)
  p <- ggplot(data=temp, aes_string(x="TSNE1",y="TSNE2",color=LMO)) + geom_point() +
    scale_color_gradient(low="lightgrey",high="red") +
    theme_void()
  print(p)
  dev.off()
}
```

```{r}
bar_filtered <- P14_1_bar_table[,1:6]
#bar_filtered <- bar_filtered[,-bad_barcodes]
```

```{r}
bar.table_sweep.list <- list()
n <- 0
for (q in seq(0.01, 0.99, by=0.02)) {
  print(q)
  n <- n + 1
  bar.table_sweep.list[[n]] <- classifyCells(bar_filtered, q=q)
  names(bar.table_sweep.list)[n] <- paste("q=",q,sep="")
}
```

```{r}
threshold.results1 <- findThresh(call.list=bar.table_sweep.list)
ggplot(data=threshold.results1$res, aes(x=q, y=Proportion, color=Subset)) + geom_line() + theme(legend.position = "none") +
  geom_vline(xintercept=threshold.results1$extrema, lty=2) + scale_color_manual(values=c("red","black","blue"))
```

```{r}
round1.calls <- classifyCells(bar_filtered, 
                              q=findQ(threshold.results1$res, 
                                      threshold.results1$extrema))
```

View classifications
```{r}
tsne$Classification <- "Singlet"
tsne$Classification[which(round1.calls[rownames(tsne)]=="Doublet")] <- "Doublet"
tsne$Classification[which(round1.calls[rownames(tsne)]=="Negative")] <- "Negative"
```

```{r}
table(tsne$Classification)
```


```{r}
tsne_classification <- ggplot(tsne, aes(x=TSNE1, y=TSNE2)) +
  geom_point(size=0.25, aes(color=Classification)) +
  theme_void()
ggsave(tsne_classification, file="../figs/P14_1_classification_tsne.png",
       dpi=600, width=8, height=5.5)
tsne_classification
```

Reclassify negative cells
```{r}
bar_full <- P14_1_bar_table[,1:6]
reclass.cells <- findReclassCells(bar_full,
                                  unique(names(round1.calls)[which(round1.calls=="Negative")]))
reclass.res <- rescueCells(bar_full, round1.calls[rownames(bar_full)], reclass.cells)
```

```{r}
ggplot(reclass.res[-1, ], aes(x=ClassStability, y=MatchRate_mean)) + 
    geom_point() + xlim(c(nrow(reclass.res)-1,1)) + 
    ylim(c(0,1.05)) +
    geom_errorbar(aes(ymin=MatchRate_mean-MatchRate_sd, ymax=MatchRate_mean+MatchRate_sd), width=.1) +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1], color="red") +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1]+3*reclass.res$MatchRate_sd[1], color="red",lty=2) +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1]-3*reclass.res$MatchRate_sd[1], color="red",lty=2)
```

```{r}
final.calls <- round1.calls
rescue.ind <- which(reclass.cells$ClassStability >= 4) ## Note: Value will be dataset-specific
final.calls[rownames(reclass.cells)[rescue.ind]] <- reclass.cells$Reclassification[rescue.ind]
```

Recheck tSNE
```{r}
tsne$Classification <- "Singlet"
tsne$Classification[which(final.calls[rownames(tsne)]=="Doublet")] <- "Doublet"
tsne$Classification[which(final.calls[rownames(tsne)]=="Negative")] <- "Negative"
```

```{r}
table(tsne$Classification)
```

```{r}
tsne_classification <- ggplot(tsne, aes(x=TSNE1, y=TSNE2)) +
  geom_point(size=0.25, aes(color=Classification)) +
  theme_void()
ggsave(tsne_classification, file="../figs/P14_1_classification_tsne.png",
       dpi=600, width=8, height=5.5)
tsne_classification
```

Plot annotations on tsne plot
```{r}
samples <- unique(final.calls)
plotSampleTSNE <- function(sample){
  data <- tsne
  data$Sample <- "Other"
  data$Sample[which(final.calls[rownames(tsne)]==sample)] <- sample
  sample_plot <- ggplot(data, aes(x=TSNE1, y=TSNE2)) +
    geom_point(size=0.25, alpha=0.5, aes(color=Sample)) +
    scale_color_manual(values=c("red", "lightgrey")) +
    theme_void()
  ggsave(sample_plot, file=paste0("../figs/LMOs/P14_1/Classifications/", 
         sample,".png"), width=5, height=3.2)
}
```

```{r}
lapply(samples, plotSampleTSNE)
```

```{r}
P14_1_calls <- final.calls
```

## P14_2
```{r}
tsne <- barTSNE(P14_2_bar_table[,1:6])
```

```{r}
plot(tsne$TSNE1, tsne$TSNE2, pch=20)
```

```{r}
temp <- P14_2_bar_table
temp$TSNE1 <- tsne$TSNE1
temp$TSNE2 <- tsne$TSNE2
write.csv(temp, file="../output/P14_2_barcode_tsne.csv", quote=F)
```

```{r}
temp <- as.matrix(tsne[3:8])
temp[temp<0] <- 0
temp <- as.data.frame(temp)
temp$TSNE1 <- tsne$TSNE1
temp$TSNE2 <- tsne$TSNE2


for (LMO in colnames(tsne)[3:8]){
  png(filename = paste0("../figs/LMOs/P14_2/",LMO,".png"), width = 700, height = 600)
  p <- ggplot(data=temp, aes_string(x="TSNE1",y="TSNE2",color=LMO)) + geom_point() +
    scale_color_gradient(low="lightgrey",high="red") +
    theme_void()
  print(p)
  dev.off()
}
```

```{r}
bar_filtered <- P14_2_bar_table[,1:6]
#bar_filtered <- bar_filtered[,-bad_barcodes]
```

```{r}
bar.table_sweep.list <- list()
n <- 0
for (q in seq(0.01, 0.99, by=0.02)) {
  print(q)
  n <- n + 1
  bar.table_sweep.list[[n]] <- classifyCells(bar_filtered, q=q)
  names(bar.table_sweep.list)[n] <- paste("q=",q,sep="")
}
```

```{r}
threshold.results1 <- findThresh(call.list=bar.table_sweep.list)
ggplot(data=threshold.results1$res, aes(x=q, y=Proportion, color=Subset)) + geom_line() + theme(legend.position = "none") +
  geom_vline(xintercept=threshold.results1$extrema, lty=2) + scale_color_manual(values=c("red","black","blue"))
```

```{r}
round1.calls <- classifyCells(bar_filtered, 
                              q=findQ(threshold.results1$res, 
                                      threshold.results1$extrema))
```

View classifications
```{r}
tsne$Classification <- "Singlet"
tsne$Classification[which(round1.calls[rownames(tsne)]=="Doublet")] <- "Doublet"
tsne$Classification[which(round1.calls[rownames(tsne)]=="Negative")] <- "Negative"
```

```{r}
table(tsne$Classification)
```


```{r}
tsne_classification <- ggplot(tsne, aes(x=TSNE1, y=TSNE2)) +
  geom_point(size=0.25, aes(color=Classification)) +
  theme_void()
ggsave(tsne_classification, file="../figs/P14_2_classification_tsne.png",
       dpi=600, width=8, height=5.5)
tsne_classification
```

Reclassify negative cells
```{r}
bar_full <- P14_2_bar_table[,1:6]
reclass.cells <- findReclassCells(bar_full,
                                  unique(names(round1.calls)[which(round1.calls=="Negative")]))
reclass.res <- rescueCells(bar_full, round1.calls[rownames(bar_full)], reclass.cells)
```

```{r}
ggplot(reclass.res[-1, ], aes(x=ClassStability, y=MatchRate_mean)) + 
    geom_point() + xlim(c(nrow(reclass.res)-1,1)) + 
    ylim(c(0,1.05)) +
    geom_errorbar(aes(ymin=MatchRate_mean-MatchRate_sd, ymax=MatchRate_mean+MatchRate_sd), width=.1) +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1], color="red") +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1]+3*reclass.res$MatchRate_sd[1], color="red",lty=2) +
    geom_hline(yintercept = reclass.res$MatchRate_mean[1]-3*reclass.res$MatchRate_sd[1], color="red",lty=2)
```

```{r}
final.calls <- round1.calls
rescue.ind <- which(reclass.cells$ClassStability >= 11) ## Note: Value will be dataset-specific
final.calls[rownames(reclass.cells)[rescue.ind]] <- reclass.cells$Reclassification[rescue.ind]
```

Recheck tSNE
```{r}
tsne$Classification <- "Singlet"
tsne$Classification[which(final.calls[rownames(tsne)]=="Doublet")] <- "Doublet"
tsne$Classification[which(final.calls[rownames(tsne)]=="Negative")] <- "Negative"
```

```{r}
table(tsne$Classification)
```

```{r}
tsne_classification <- ggplot(tsne, aes(x=TSNE1, y=TSNE2)) +
  geom_point(size=0.25, aes(color=Classification)) +
  theme_void()
ggsave(tsne_classification, file="../figs/P14_2_classification_tsne.png",
       dpi=600, width=8, height=5.5)
tsne_classification
```

Plot annotations on tsne plot
```{r}
samples <- unique(final.calls)
plotSampleTSNE <- function(sample){
  data <- tsne
  data$Sample <- "Other"
  data$Sample[which(final.calls[rownames(tsne)]==sample)] <- sample
  sample_plot <- ggplot(data, aes(x=TSNE1, y=TSNE2)) +
    geom_point(size=0.25, alpha=0.5, aes(color=Sample)) +
    scale_color_manual(values=c("red", "lightgrey")) +
    theme_void()
  ggsave(sample_plot, file=paste0("../figs/LMOs/P14_2/Classifications/", 
         sample,".png"), width=5, height=3.2)
}
```

```{r}
lapply(samples, plotSampleTSNE)
```

```{r}
P14_2_calls <- final.calls
```

# Put annotations into seurat
## P3_P7_1
```{r}
P3_P7_1_seurat$Barcode <- P3_P7_1_calls

P3_P7_1_seurat$Age <- "NA"
P3_P7_1_seurat$Age[P3_P7_1_seurat$Barcode %in% c("Bar1", "Bar2", "Bar3",
                                                 "Bar4", "Bar5", "Bar6")] <- "P3"
P3_P7_1_seurat$Age[P3_P7_1_seurat$Barcode %in% c("Bar7", "Bar8", "Bar9",
                                                 "Bar10", "Bar11", "Bar12")] <- "P7"

P3_P7_1_seurat$Oxygen <- "NA"
P3_P7_1_seurat$Oxygen[P3_P7_1_seurat$Barcode %in% c("Bar1", "Bar2", "Bar3",
                                                    "Bar7", "Bar8", "Bar9")] <- "Normoxia"
P3_P7_1_seurat$Oxygen[P3_P7_1_seurat$Barcode %in% c("Bar4", "Bar5", "Bar6",
                                                    "Bar10", "Bar11", "Bar12")] <- "Hyperoxia"
```

## P3_P7_2
```{r}
P3_P7_2_seurat$Barcode <- P3_P7_2_calls

P3_P7_2_seurat$Age <- "NA"
P3_P7_2_seurat$Age[P3_P7_2_seurat$Barcode %in% c("Bar1", "Bar2", "Bar3",
                                                 "Bar4", "Bar5", "Bar6")] <- "P3"
P3_P7_2_seurat$Age[P3_P7_2_seurat$Barcode %in% c("Bar7", "Bar8", "Bar9",
                                                 "Bar10", "Bar11", "Bar12")] <- "P7"

P3_P7_2_seurat$Oxygen <- "NA"
P3_P7_2_seurat$Oxygen[P3_P7_2_seurat$Barcode %in% c("Bar1", "Bar2", "Bar3",
                                                    "Bar7", "Bar8", "Bar9")] <- "Normoxia"
P3_P7_2_seurat$Oxygen[P3_P7_2_seurat$Barcode %in% c("Bar4", "Bar5", "Bar6",
                                                    "Bar10", "Bar11", "Bar12")] <- "Hyperoxia"
```

## P3_2
```{r}
P3_2_seurat$Barcode <- P3_2_calls

P3_2_seurat$Age <- "NA"
P3_2_seurat$Age[P3_2_seurat$Barcode %in% c("Bar1", "Bar2", "Bar3",
                                           "Bar4", "Bar5", "Bar6")] <- "P3"

P3_2_seurat$Oxygen <- "NA"
P3_2_seurat$Oxygen[P3_2_seurat$Barcode %in% c("Bar1", "Bar2", "Bar3")] <- "Normoxia"
P3_2_seurat$Oxygen[P3_2_seurat$Barcode %in% c("Bar4", "Bar5", "Bar6")] <- "Hyperoxia"
```

## P7_1
```{r}
P7_1_seurat$Barcode <- P7_1_calls

P7_1_seurat$Age <- "NA"
P7_1_seurat$Age[P7_1_seurat$Barcode %in% c("Bar1", "Bar2", "Bar3",
                                           "Bar4", "Bar5", "Bar6")] <- "P7"

P7_1_seurat$Oxygen <- "NA"
P7_1_seurat$Oxygen[P7_1_seurat$Barcode %in% c("Bar1", "Bar2", "Bar3")] <- "Normoxia"
P7_1_seurat$Oxygen[P7_1_seurat$Barcode %in% c("Bar4", "Bar5", "Bar6")] <- "Hyperoxia"
```

## P14_1
```{r}
P14_1_seurat$Barcode <- P14_1_calls

P14_1_seurat$Age <- "NA"
P14_1_seurat$Age[P14_1_seurat$Barcode %in% c("Bar1", "Bar2", "Bar3",
                                             "Bar4", "Bar5", "Bar6")] <- "P14"

P14_1_seurat$Oxygen <- "NA"
P14_1_seurat$Oxygen[P14_1_seurat$Barcode %in% c("Bar1", "Bar2", "Bar3")] <- "Normoxia"
P14_1_seurat$Oxygen[P14_1_seurat$Barcode %in% c("Bar4", "Bar5", "Bar6")] <- "Hyperoxia"
```

## P14_2
```{r}
P14_2_seurat$Barcode <- P14_2_calls

P14_2_seurat$Age <- "NA"
P14_2_seurat$Age[P14_2_seurat$Barcode %in% c("Bar1", "Bar2", "Bar3",
                                             "Bar4", "Bar5", "Bar6")] <- "P14"

P14_2_seurat$Oxygen <- "NA"
P14_2_seurat$Oxygen[P14_2_seurat$Barcode %in% c("Bar1", "Bar2", "Bar3")] <- "Normoxia"
P14_2_seurat$Oxygen[P14_2_seurat$Barcode %in% c("Bar4", "Bar5", "Bar6")] <- "Hyperoxia"
```

# Merge seurat objects
```{r}
seurat <- merge(P3_P7_1_seurat, c(P3_P7_2_seurat, P3_2_seurat,
                                  P7_1_seurat, P14_1_seurat,
                                  P14_2_seurat),
                add.cell.ids=c("P3_P7_1", "P3_P7_2", "P3_2",
                               "P7_1", "P14_1", "P14_2"))
```

```{r}
saveRDS(seurat, file="../output/merged_data_unfiltered.rds")
```