covid-sc.Rmd

---
title: "Corona Time"
output: html_notebook
---

```{r}
library(Seurat)
library(dplyr)
library(Matrix)

set.seed(6112)

```


Load a sample for experimenting and figuring things out
```{r}
rdata <- readRDS(file="covdata/GSM4557327_555_1_cell.counts.matrices.rds", refhook = NULL)
```



```{r}

datalist = list()

index <- 1

key <- list()

gnames <- list()

sampnames <- character(13)




for (n in list.files("covdata")) {
  
  
  patient <- readRDS(file=paste("covdata/", n, sep=""), refhook = NULL)
  
  patient <- patient$exon
  
  sample <- substr(n, 1, 10)
  
  sampnames[index] <- sample
  
  
  key[[letters[index]]] <- sample
  
  
  #cnamesedit <- character(length(patient$exon@Dimnames[[2]]))
  
  
  # for (i in 1:length(patient$exon@Dimnames[[2]])) {
  #   cnamesedit[i] <- paste(letters[index], patient$exon@Dimnames[[2]][i], sep = "")
  # }
  
  
  patient@Dimnames[[2]] <- gsub("(^[ATCG]{12}$)", paste0(letters[index], "\\1"), patient@Dimnames[[2]])
  
  
  
  gnames[[index]] <- patient@Dimnames[[1]]
  
  
  
  datalist[[sample]] <- patient
  
  index <- index + 1
}
```



```{r}


allgnames <- Reduce(intersect, gnames)
sindex <- 1

for (sample in datalist) {
  indices <- match(allgnames, sample@Dimnames[[1]])
  
  datalist[[sampnames[sindex]]] <- sample[indices,]
  
  #e <- sample[indices,]
  
  sindex <- sindex + 1
}

#gnames <- datalist$GSM4557327@Dimnames[[1]]

#mnames <- datalist$GSM4557327@Dimnames[[2]]


```


```{r}

alldata <- do.call(cbind, datalist)


###
saveRDS(alldata, file = "alldata.rds")
###

```



```{r}
# exongenes <- rdata$exon@Dimnames[[1]]
# introngenes <- rdata$intron@Dimnames[[1]]
# spanninggenes <- rdata$spanning@Dimnames[[1]]
# 
# 
# 
# exons <- datalist$GSM4557327@x




```


Calculate different aspects of data to help determine data cut offs
```{r}
counts_per_cell <- Matrix::colSums(alldata)
counts_per_gene <- Matrix::rowSums(alldata)

genes_per_cell <- Matrix::colSums(alldata>0)
cells_per_gene <- Matrix::rowSums(alldata>0)

gene_count_ratio <- genes_per_cell / counts_per_cell


edited_cell_counts <- counts_per_cell[-which(gene_count_ratio > 0.92)]
edited_cell_genes <- genes_per_cell[-which(gene_count_ratio > 0.92)]

edited_gene_count_ratio <- edited_cell_genes / edited_cell_counts

editdata <- alldata[,-which(gene_count_ratio > 0.92)]

```


```{r}
geneplot <- plot(sort(genes_per_cell), xlab='cell', log='y', main='genes per cell (ordered)')

geneplotedit <- plot(sort(edited_cell_genes), xlab='cell', log='y', main='genes per cell (ordered)')

abline(h = 100, col = "red")
abline(h = 3000, col = "red")
abline(h = 1750, col = "red")


genecellplot <- plot(sort(cells_per_gene), xlab='gene', log='y', main='cells per gene (ordered)')


```

```{r}
plot(gene_count_ratio)
plot(edited_gene_count_ratio)
```



```{r}

#counts_per_cell <- Matrix::colSums(alldata)

###genes_per_cell <- Matrix::colSums(alldata[,16833:16834]>0)  
###16833 is a point on the strange line of points along line y = x

#genes_per_cell <- Matrix::colSums(alldata>0)

###plot(counts_per_cell, genes_per_cell, col=as.factor(c(rep(1, length(counts_per_cell)), rep(2,))), cex=1, xlim = c(0,20000), ylim = c(0, 5000))

plot(counts_per_cell, genes_per_cell, col=as.factor(c(rep(1, length(counts_per_cell)), rep(2,))), cex=1, log = "xy")


abline(a = 0, b = 1, col = "blue")

```


```{r}
plot(edited_cell_counts, edited_cell_genes, log="xy")
```



```{r}
hist(log10(counts_per_cell+1),main='counts per cell',col='wheat', breaks = 50)
hist(log10(genes_per_cell+1), main='genes per cell', col='wheat', breaks = 50)
plot(counts_per_cell, genes_per_cell, log='xy', col=as.factor(substring(names(counts_per_cell), 1,1)), cex=1)
#plot(counts_per_cell, genes_per_cell, col=as.factor(c(rep(1, length(counts_per_cell)), rep(2,))), cex=1)
title('counts vs genes per cell')


hist(log10(edited_cell_counts+1),main='counts per cell',col='wheat', breaks = 50)
hist(log10(edited_cell_genes+1), main='genes per cell', col='wheat', breaks = 50)
plot(edited_cell_counts, edited_cell_genes, log='xy', col=as.factor(substring(names(edited_cell_counts), 1,1)), cex=1)

title('counts vs genes per cell')
#plot(cells_per_gene, counts_per_gene, log="xy", cex=1)
```

Create patient metadata dataframe
```{r}

status <- gsub("^[a-g]{1}[ATCG]{12}$", "covid", editdata@Dimnames[[2]])

status <- gsub("^[h-m]{1}[ATCG]{12}$", "healthy", status)


meta <- data.frame("status" = status, row.names = editdata@Dimnames[[2]])

```


Create Seurat object
```{r}
seur <- CreateSeuratObject(counts = editdata, project = "covidsc", min.cells = 5, min.features = 100, meta.data = meta)

#max ncount 12500
#min ncount 125
#max nfeature 2500
#min nfeature 200

#nfeatures = number of genes expressed by each cell
#ncounts = number of total molecules in a cell
```

```{r}
seur[["percent.mt"]] <- PercentageFeatureSet(seur, pattern = "^MT-")

VlnPlot(seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
```


```{r}
plot1 <- FeatureScatter(seur, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

plot1
plot2

```




Abnormally large amounts of mitochondria likely means cell is dying (mitochondrial contamination)
```{r}
mtplot <- plot(sort(seur$percent.mt), xlab='% mito', log='y', main='mito percent per cell (ordered)')

abline(h=10, col = "red")
```


Apply cut offs to reduce outliers
```{r}
seur <- subset(seur, subset = 
                 nCount_RNA > 1000 
               & nCount_RNA < 12500 
               & nFeature_RNA > 150 
               & nFeature_RNA < 2500 
               & percent.mt < 10)


#seur <- subset(seur, subset = nCount_RNA > 125 & nCount_RNA < 12500 & percent.mt < 5)

#16199 
#ncount min - 15916 | dif 283
#ncount max - 15904 | dif 12
#nfeature min - 10082 | dif 5822
#nfeature max - 9996 | dif 86
#mito - 8711 | dif 1287

```


```{r}
#alldata - 150370 cells
#editdata - 147116 cells
#seur initial - 136918 cells
#seur after subset - 39512 cells
```


```{r}
plot3 <- FeatureScatter(seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

plot3
```




Normalize the data with default methods
```{r}
seur <- NormalizeData(seur, normalization.method = "LogNormalize", scale.factor = 10000)
```








```{r}
seur <- FindVariableFeatures(seur, selection.method = "vst", nfeatures = 2000)

length(VariableFeatures(seur))


top10 <- head(VariableFeatures(seur), 10)

top10


varplot <- VariableFeaturePlot(seur)
varplot2 <- LabelPoints(plot = varplot, points = top10, repel = TRUE)

varplot
varplot2

```



scale data
```{r}

#todo
#understand what is going on
#research and figure out good high low cut offs for noise
#fix code below and fine tune to fit covid data

allgenes <- rownames(seur)
seur <- ScaleData(seur, features = allgenes)
```




#####################################


PCA
```{r}
seur <- RunPCA(seur, features = VariableFeatures(object = seur))
```

```{r}
print(seur[["pca"]], dims = 1:5, nfeatures = 5)
```

```{r}
VizDimLoadings(seur, dims = 1:2, reduction = "pca")
```

```{r}
DimPlot(seur, reduction = "pca")
```

```{r}
DimHeatmap(seur, dims = 1, cells = 500, balanced = TRUE)
```

```{r}
DimHeatmap(seur, dims = 1:9, cells = 500, balanced = TRUE)
```

```{r}
ElbowPlot(seur, ndims = 50)

#PC cutoff around 15 or 22
```

Find neighbors
```{r}
seur <- FindNeighbors(seur, reduction = "pca", dims = 1:22)
```



Default Louvain Clustering
```{r}
seur <- FindClusters(seur, resolution = 0.4)



#0.4 - 0.9303
#0.5 - 0.9201
#0.6 - 0.9116
#0.7 - 0.9038
#0.8 - 0.8960 -- umpa is good
#0.9 - 0.8894
#1.2 - 0.8711


head(Idents(seur), 5)


seur <- RunUMAP(seur, dims = 1:22)

DimPlot(seur, reduction = "umap")
```


SLM clustering
```{r}
seur <- FindClusters(seur, resolution = 0.8, algorithm = 3)

head(Idents(seur), 5)


seur <- RunUMAP(seur, dims = 1:22)

DimPlot(seur, reduction = "umap")
```



Louvain
```{r}
seur <- RunUMAP(seur, dims = 1:22)

DimPlot(seur, reduction = "umap")
```


SLM
```{r}
seur <- RunUMAP(seur, dims = 1:22)

DimPlot(seur, reduction = "umap")
```