library(Seurat)
library(dplyr)

tam.data <- Read10X(data.dir="TAMs/patient-1/")
tam <- CreateSeuratObject(counts=tam.data, project="TAMs", min.cells=10, min.features=800)
tam

tam[["percent.mt"]] <- PercentageFeatureSet(tam, pattern="^mt-")
plot <- FeatureScatter(tam, feature1="nFeature_RNA", feature2="percent.mt")
plot
plot <- FeatureScatter(tam, feature1="nCount_RNA", feature2="percent.mt")
plot

tam[["reliable"]] <- tam[["nCount_RNA"]]<=65000 & tam[["percent.mt"]]<=12
tam <- subset(tam, subset=reliable)
tam <- NormalizeData(tam)

all.genes <- rownames(tam)
tam <- ScaleData(tam, features=all.genes)
tam <- FindVariableFeatures(tam, nfeatures=2000)
ccycle <- read.csv("TAMs/GO-mouse-cell-cycle.txt", sep="\t", stringsAsFactors=FALSE, check.names = FALSE)
VariableFeatures(tam) <- setdiff(VariableFeatures(tam), ccycle$Symbol)
tam <- RunPCA(tam, features=VariableFeatures(tam), ndims.print=1:3, nfeatures.print=10)
tam <- FindNeighbors(tam, dims=1:35)
tam <- FindClusters(tam, resolution=1.5) # resolution controls the size of the communities
tam <- RunUMAP(tam, dims=1:35)
DimPlot(tam, reduction="umap")
tam <- RunTSNE(tam, dims=1:35)
DimPlot(tam, reductio="tsne")
DimPlot(tam, reduction="pca")

tam.markers <- FindAllMarkers(tam, only.pos=TRUE, min.pct=0.25, logfc.threshold=0.25)
tam.markers %>%
  group_by(cluster) %>%
  slice_max(n=2, order_by=avg_log2FC)

top10 <- tam.markers %>%
  group_by(cluster) %>%
  slice_max(n=10, order_by=avg_log2FC)

DoHeatmap(tam, features=top10$gene) + NoLegend()

VlnPlot(tam, features=c("Tmem86a","Mki67","Ccr2","Cxcr4","Mrc1","Tnip3","Ccl7"))


library(destiny)
pca.proj <- Embeddings(tam, reduction="pca")
dm <- DiffusionMap(pca.proj[, 1:35])
plot(dm)
cluster.num <- as.numeric(Idents(tam))
cols <- rainbow(max(cluster.num))
plot(x=dm@eigenvectors[,1], y=dm@eigenvectors[,2], col=cols[cluster.num], pch=20)
legend(x="topleft",legend=1:max(cluster.num),pch=20,col=cols[1:max(cluster.num)])

library(slingshot)
sl <- slingshot(dm@eigenvectors[,1:15], clusterLabels=cluster.num, start.clus=1)
sds <- as.SlingshotDataSet(sl)
lines(sds, type="c", lwd=2)
lin <- getLineages(dm@eigenvectors[,1:15], clusterLabels=cluster.num, start.clus=1)
lines(SlingshotDataSet(lin), lwd = 3, col = 'gray')