In MyersGroup/testisAtlas: testisAtlas: R functions and analysis scripts used in the analysis of testis single cell RNAseq

print_tsne("Prdm9", curve = T) + theme(panel.grid.major = element_blank(),
                                       panel.grid.minor = element_blank(),
                                       axis.line = element_blank(),
                                       axis.title.x=element_blank(),
                                       axis.title.y=element_blank(),
                                       axis.text.x=element_blank(),
                                       axis.text.y=element_blank())

#devtools::install_github("thomasp85/tweenr")
#devtools::install_github("dgrtwo/gganimate")

library(tweenr)
library(gganimate)
library(ggplot2)
library(viridis)

devtools::load_all() # aka library(testisAtlas)

load2("../data/cache")

tmp <- merge(cell_data[somatic4==FALSE,
                   c("cell","PseudoTime","PseudoTime1","PseudoTime2","Tsne1_QC1", "Tsne2_QC1","group"),
                   with=FALSE],
               sda_predict("Prdm9", name_extension = ""))

# define two key arrangements of points

ts <- list(data.frame(tmp[,.("x"=(Tsne1_QC1-mean(Tsne1_QC1))/sd(Tsne1_QC1),
                             "y"=(Tsne2_QC1-mean(Tsne2_QC1))/sd(Tsne2_QC1), pt=Prdm9)]),
           # comment out to collapse points to pseudotime line
           #data.frame(tmp[,.("x"=(PseudoTime1-mean(Tsne1_QC1))/sd(Tsne1_QC1),
           #                 "y"=(PseudoTime2-mean(Tsne2_QC1))/sd(Tsne2_QC1), pt=Prdm9)]),
           data.frame(tmp[,.("x"=scale(-PseudoTime), "y"=Prdm9/(4*sd(Prdm9))-1.5, pt=Prdm9)]))


tf <- tween_states(ts, tweenlength = 10, statelength = 1, 
                   ease = c('cubic-in-out'), 
                   nframes = 99)


# pseudotime curve data
curve_data <- data.frame(x = (principal_curves[["df_9"]]$s[principal_curves[["df_9"]]$tag, 1]
                              - mean(tmp$Tsne1_QC1)) / sd(tmp$Tsne1_QC1),
                         y = (principal_curves[["df_9"]]$s[principal_curves[["df_9"]]$tag, 2]
                              - mean(tmp$Tsne2_QC1)) / sd(tmp$Tsne2_QC1))

# simplify the curve: use every 100th point (-> 168 points)
curve_data <- curve_data[seq(from = 1,to = 16715, by = 100),]

# check still looks ok
# ggplot(curve_data, aes(x, y)) + geom_path()

# show for 10 frames
curve_data <- do.call("rbind", replicate(10, curve_data, simplify = FALSE))

curve_data$.frame <- rep(1:10, each=168)
curve_data$.fade <- rep(seq(1,0.1,-0.1), each=168)


# Animate with gganimate
p <- ggplot(data=tf, aes(x=x, y=y, colour=pt)) + 
  geom_point(aes(frame = .frame), size=0.5) + 
  scale_colour_viridis(direction = -1) +
  theme_minimal(base_size = 21) +
  labs(colour="Prdm9 Expression") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.x=element_blank(),
        axis.text.y=element_blank()) +
  theme(legend.position="bottom")

# add curve path to main plot
p <- p + geom_path(data = curve_data,
                   aes(frame = .frame, alpha=.fade),
                   size = 0.5,
                   colour = "Black",
                   arrow = arrow(angle = 15, ends = "first", type = "closed")
                   ) +
          guides(alpha=FALSE)

animation::ani.options(interval = 1/15)

gganimate(p, "../results/tSNE_pseudotime.gif", title_frame = F, ani.width = 1000, ani.height = 1000)

ts <- list(data.frame(melt_genes(c("Prdm9","Gapdhs"))),
           data.frame(melt_genes(c("Prdm9","Gapdhs"), predict = T)),
           data.frame(melt_genes(c("Prdm9","Gapdhs"))))

ts[[1]]$Type <- "Raw Normalised"
ts[[2]]$Type <- "Imputed"
ts[[3]]$Type <- "Raw Normalised"

ts[[1]]$col <- RColorBrewer::brewer.pal(3, "Set1")[1]
ts[[2]]$col <- RColorBrewer::brewer.pal(3, "Set1")[2]
ts[[3]]$col <- RColorBrewer::brewer.pal(3, "Set1")[1]

tf <- tween_states(ts, tweenlength = 10, statelength = 1, 
                   ease = c('cubic-in-out'), 
                   nframes = 99)

p <- ggplot(data=tf, aes(x=-PseudoTime, y=Expression)) + 
  geom_point(aes(frame = .frame, colour=col), size=0.5) + 
  scale_colour_identity() +
  theme_minimal(base_size = 25) +
  theme(legend.position = "bottom") +
  facet_wrap(~Gene, ncol = 1, scales = "free_y") +
  ggtitle("Raw Normalised (red) vs Imputed (blue)") +
  xlab("PseudoTime")

# p <- ggplot(data=tf, aes(x=-PseudoTime, y=Gapdhs)) + 
#   geom_point(aes(frame = .frame), size=0.5) + 
#   theme_minimal()

gganimate(p, "../results/Impuation.gif", title_frame = F, ani.width = 1400, ani.height = 1000)

MyersGroup/testisAtlas documentation built on Nov. 29, 2020, 9:21 p.m.