***
library(Seurat)
library(dplyr)

# FILTERED matrix from CellRanger
all.data <- Read10X(data.dir = "path")
all <- CreateSeuratObject(raw.data = all.data, min.cells = 0, min.genes = 0, project = "10X_all")

# Filtering - first find the mitochondria genes and stash info in MetaData
mito.genes <- grep(pattern = "^mt-", x = rownames(x = all@data), value = TRUE)
percent.mito <- Matrix::colSums(all@raw.data[mito.genes, ]) / Matrix::colSums(all@raw.data)
all <- AddMetaData(object = all, metadata = percent.mito, col.name = "percent.mito")
all <- FilterCells(object = all, subset.names = c("nGene", "percent.mito", "nUMI"), low.thresholds = c(200, -Inf, -Inf), high.thresholds = c(Inf, 0.25, 100000))

# Normalize the data
all <- NormalizeData(object = all, normalization.method = "LogNormalize", scale.factor = 10000)

# Find Variable Genes
all <- FindVariableGenes(object = all, mean.function = ExpMean, dispersion.function = LogVMR, x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5)
length(x = all@var.genes) ## check the length to see if reasonable

# Regress out information that shouldn't be used to cluster the cells 
all <- ScaleData(object = all, vars.to.regress = c("nUMI", "percent.mito"), do.par = TRUE, num.cores = 4)

# Identify PCs 
all <- RunPCA(object = all, pc.genes = all@var.genes, do.print = TRUE, pcs.print = 1:5, genes.print = 5, pcs.compute = 200)
all <- ProjectPCA(object = all, do.print = FALSE, pcs.store = 200)

# JackStraw to help identify the number of PCs to use
all <- JackStraw(object = all, num.pc = 200, num.replicate = 25, display.progress = TRUE,  do.par = TRUE, num.cores = 4)
jackStraw <- JackStrawPlot(object = all, PCs = 1:200)
ggsave(paste("jackStraw.pdf", sep=""), width=15, height=100, jackStraw, device = "jpeg", limitsize = FALSE)

# Cluster the cells based on the number of PCs selected 
all <- FindClusters(object = all, reduction.type = "pca", dims.use = 1:179, resolution = 18.0, print.output = FALSE, save.SNN = TRUE, force.recalc = TRUE)

# tSNE for visualization
all <- RunTSNE(object = all, dims.use = 1:179)

# Find all markers
all.markers <- FindAllMarkers(object = all, only.pos = TRUE, min.pct = 0.25, logfc.threshold = log(0.5))
write.csv(all.markers, file="allMarkers_GRCz11_179_18.csv")
***

# UMAP and metadata to matrix

all <- RunUMAP(object = all, reduction.use = "pca", dims.use = 1:179, min_dist = 0.555, n_neighbors = 37)


```{r atlas plots, fig.height=20, fig.width=20}
# TSNE commands here to allow for an alternative view
TSNEPlot(object = all, do.label = TRUE, no.legend = T)

# to plot with UMAP coordinates the below commands are needed
DimPlot(object = all, reduction.use = "umap", do.label = T, no.legend = T, do.return = F, vector.friendly = F, pt.size = 0.1, label.size = 6)

```