== Making diffusion maps with destiny ==

Seurat v3 doesn't have a function to make diffusion maps anymore. The alternatives are to 1) run the analysis in python using scanpy or 2) make a single cell object, i.e. from the Seurat object, and run the "DiffusionMap" function from the destiny package. Example commands for the second option are shown below. 
 
  * **Making a single cell object from a Seurat object** 
{{{
library(Seurat) 
library(SingleCellExperiment)
sce <- as.SingleCellExperiment(seuratObject)
#this has the cell classification
table(sce$ident)
}}}

  * **Running diffusion map** 
{{{
#this step may take a long time (days) or not finish. It is recommend to send it to the cluster as a script that reads the Seurat or the single cell object, runs DiffusionMap, and saves the object. 
library(destiny)
dm <- DiffusionMap(sce, verbose = TRUE)
}}}

  * **Plotting the diffusion map** 
{{{
library(ggplot2)
cellLabels <- sce$ident
tmp <- data.frame(DC1 = eigenvectors(dm)[, 1],
                  DC2 = eigenvectors(dm)[, 2],
                  DC3 = eigenvectors(dm)[, 3],
                  DC4 = eigenvectors(dm)[, 4],
                  Samples = cellLabels)
pdf("./DC1_DC2.pdf", w=11, h=8.5)
ggplot(tmp, aes(x = DC1, y = DC2, colour = Samples)) +
  geom_point()  + 
  xlab("Diffusion component 1") + 
  ylab("Diffusion component 2") +
  theme_classic()
dev.off()
}}}
  * **Plotting cell progression along the diffusion map components** 
{{{
sce$pseud_dm1 <- rank(eigenvectors(dm)[,1])      # rank cells by their dpt dm1
sce$pseud_dm2 <- rank(eigenvectors(dm)[,2])      # rank cells by their dpt dm2
sce$pseud_dm1R <- rank(-eigenvectors(dm)[,1])    # rank cells by their dpt dm1 reverse order
sce$pseud_dm2R <- rank(-eigenvectors(dm)[,2])    # rank cells by their dpt dm2 reverse order

SortedDM1 <- data.frame(DM1Sort = as.data.frame(colData(sce))$pseud_dm1,
                        Samples = as.data.frame(colData(sce))$ident)
SortedDM2 <- data.frame(DM2Sort = as.data.frame(colData(sce))$pseud_dm2,
                        Samples = as.data.frame(colData(sce))$ident)
SortedDM1R <- data.frame(DM1SortR = as.data.frame(colData(sce))$pseud_dm1R,
                         Samples = as.data.frame(colData(sce))$ident)
SortedDM2R <- data.frame(DM2SortR = as.data.frame(colData(sce))$pseud_dm2R,
                         Samples = as.data.frame(colData(sce))$ident)

ggplot(SortedDM1, aes(x=SortedDM1[,1], y=Samples,color=Samples)) +
  geom_jitter() + xlab("Diffusion component 1 (DC1)") + ylab("Samples") +
  ggtitle("Cells ordered by DC1")
ggplot(SortedDM2, aes(x=SortedDM2[,1], y=Samples,color=Samples)) +
  geom_jitter() + xlab("Diffusion component 2 (DC2)") + ylab("Samples") +
  ggtitle("Cells ordered by DC2")
  
ggplot(SortedDM1R, aes(x=SortedDM1R[,1], y=Samples,color=Samples)) +
  geom_jitter() + xlab("Minus Diffusion component 1 (DC1)") + ylab("Samples") +
  ggtitle("Cells ordered by reversed DC1")
ggplot(SortedDM2R, aes(x=SortedDM2R[,1], y=Samples,color=Samples)) +
  geom_jitter() + xlab("Minus Diffusion component 2 (DC2)") + ylab("Samples") +
  ggtitle("Cells ordered by reversed DC2")
}}}

  * **Make and interactive 2D and 3D diffusion map figure** 
{{{
library(plotly)

#interactive 2D
p2 = plot_ly(x=tmp$DC1, y=tmp$DC2, type="scatter", mode="markers", color=tmp$Samples, marker.size = 0.5)
htmlwidgets::saveWidget(as_widget(p2), "Interactive2D_DiffM.html", title = "Diffusion map")

#interactive 3D
p = plot_ly(x=tmp$DC1, y=tmp$DC2, z=tmp$DC3, type="scatter3d", mode="markers", color=tmp$Samples, marker = list(size = 2 ))
htmlwidgets::saveWidget(as_widget(p), "Interactive3D.html", title = "Diffusion map")

}}}

{{{

}}}