monocle的细胞轨迹图用ggplot重新做

my_plot_cell_trajectory<-function (cds, x = 1, y = 2, color_by = "State", show_tree = TRUE, 
                                   show_backbone = TRUE, backbone_color = "black", markers = NULL, 
                                   use_color_gradient = FALSE, markers_linear = FALSE, show_cell_names = FALSE, 
                                   show_state_number = FALSE, cell_size = 1.5, cell_link_size = 0.75, 
                                   cell_name_size = 2, state_number_size = 2.9, show_branch_points = TRUE, 
                                   theta = 0, ...) 
{
  requireNamespace("igraph")
  gene_short_name <- NA
  sample_name <- NA
  sample_state <- pData(cds)$State
  data_dim_1 <- NA
  data_dim_2 <- NA
  lib_info_with_pseudo <- pData(cds)
  if (is.null(cds@dim_reduce_type)) {
    stop("Error: dimensionality not yet reduced. Please call reduceDimension() before calling this function.")
  }
  if (cds@dim_reduce_type == "ICA") {
    reduced_dim_coords <- reducedDimS(cds)
  }
  else if (cds@dim_reduce_type %in% c("simplePPT", "DDRTree")) {
    reduced_dim_coords <- reducedDimK(cds)
  }
  else {
    stop("Error: unrecognized dimensionality reduction method.")
  }
  ica_space_df <- Matrix::t(reduced_dim_coords) %>% as.data.frame() %>% 
    select_(prin_graph_dim_1 = x, prin_graph_dim_2 = y) %>% 
    mutate(sample_name = rownames(.), sample_state = rownames(.))
  dp_mst <- minSpanningTree(cds)
  if (is.null(dp_mst)) {
    stop("You must first call orderCells() before using this function")
  }
  edge_df <- dp_mst %>% igraph::as_data_frame() %>% select_(source = "from", 
                                                            target = "to") %>% left_join(ica_space_df %>% select_(source = "sample_name", 
                                                                                                                  source_prin_graph_dim_1 = "prin_graph_dim_1", source_prin_graph_dim_2 = "prin_graph_dim_2"), 
                                                                                         by = "source") %>% left_join(ica_space_df %>% select_(target = "sample_name", 
                                                                                                                                               target_prin_graph_dim_1 = "prin_graph_dim_1", target_prin_graph_dim_2 = "prin_graph_dim_2"), 
                                                                                                                      by = "target")
  data_df <- t(monocle::reducedDimS(cds)) %>% as.data.frame() %>% 
    select_(data_dim_1 = x, data_dim_2 = y) %>% rownames_to_column("sample_name") %>% 
    mutate(sample_state) %>% left_join(lib_info_with_pseudo %>% 
                                         rownames_to_column("sample_name"), by = "sample_name")
  return_rotation_mat <- function(theta) {
    theta <- theta/180 * pi
    matrix(c(cos(theta), sin(theta), -sin(theta), cos(theta)), 
           nrow = 2)
  }
  rot_mat <- return_rotation_mat(theta)
  cn1 <- c("data_dim_1", "data_dim_2")
  cn2 <- c("source_prin_graph_dim_1", "source_prin_graph_dim_2")
  cn3 <- c("target_prin_graph_dim_1", "target_prin_graph_dim_2")
  data_df[, cn1] <- as.matrix(data_df[, cn1]) %*% t(rot_mat)
  edge_df[, cn2] <- as.matrix(edge_df[, cn2]) %*% t(rot_mat)
  edge_df[, cn3] <- as.matrix(edge_df[, cn3]) %*% t(rot_mat)
  markers_exprs <- NULL
  if (is.null(markers) == FALSE) {
    markers_fData <- subset(fData(cds), gene_short_name %in% 
                              markers)
    if (nrow(markers_fData) >= 1) {
      markers_exprs <- reshape2::melt(as.matrix(exprs(cds[row.names(markers_fData), 
      ])))
      colnames(markers_exprs)[1:2] <- c("feature_id", "cell_id")
      markers_exprs <- merge(markers_exprs, markers_fData, 
                             by.x = "feature_id", by.y = "row.names")
      markers_exprs$feature_label <- as.character(markers_exprs$gene_short_name)
      markers_exprs$feature_label[is.na(markers_exprs$feature_label)] <- markers_exprs$Var1
    }
  }
  if (is.null(markers_exprs) == FALSE && nrow(markers_exprs) > 
      0) {
    data_df <- merge(data_df, markers_exprs, by.x = "sample_name", 
                     by.y = "cell_id")
    if (use_color_gradient) {
      if (markers_linear) {
        g <- ggplot(data = data_df, aes(x = data_dim_1, 
                                        y = data_dim_2)) + geom_point(aes(color = value), 
                                                                      size = I(cell_size), na.rm = TRUE) + scale_color_viridis(name = paste0("value"), 
                                                                                                                               ...) + facet_wrap(~feature_label)
      }
      else {
        g <- ggplot(data = data_df, aes(x = data_dim_1, 
                                        y = data_dim_2)) + geom_point(aes(color = log10(value + 
                                                                                          0.1)), size = I(cell_size), na.rm = TRUE) + 
          scale_color_viridis(name = paste0("log10(value + 0.1)"), 
                              ...) + facet_wrap(~feature_label)
      }
    }
    else {
      if (markers_linear) {
        g <- ggplot(data = data_df, aes(x = data_dim_1, 
                                        y = data_dim_2, size = (value * 0.1))) + facet_wrap(~feature_label)
      }
      else {
        g <- ggplot(data = data_df, aes(x = data_dim_1, 
                                        y = data_dim_2, size = log10(value + 0.1))) + 
          facet_wrap(~feature_label)
      }
    }
  }
  else {
    g <- ggplot(data = data_df, aes(x = data_dim_1, y = data_dim_2))
  }
  if (show_tree) {
    g <- g + geom_segment(aes_string(x = "source_prin_graph_dim_1", 
                                     y = "source_prin_graph_dim_2", xend = "target_prin_graph_dim_1", 
                                     yend = "target_prin_graph_dim_2"), size = cell_link_size, 
                          linetype = "solid", na.rm = TRUE, data = edge_df)
  }
  if (is.null(markers_exprs) == FALSE && nrow(markers_exprs) > 
      0) {
    if (use_color_gradient) {
    }
    else {
      g <- g + geom_point(aes_string(color = color_by), 
                          na.rm = TRUE)
    }
  }
  else {
    if (use_color_gradient) {
    }
    else {
      g <- g + geom_point(aes_string(color = color_by), 
                          size = I(cell_size), na.rm = TRUE)
    }
  }
  if (show_branch_points && cds@dim_reduce_type == "DDRTree") {
    mst_branch_nodes <- cds@auxOrderingData[[cds@dim_reduce_type]]$branch_points
    branch_point_df <- ica_space_df %>% slice(match(mst_branch_nodes, 
                                                    sample_name)) %>% mutate(branch_point_idx = seq_len(n()))
    g <- g + geom_point(aes_string(x = "prin_graph_dim_1", 
                                   y = "prin_graph_dim_2"), size = 5, na.rm = TRUE, 
                        branch_point_df) + geom_text(aes_string(x = "prin_graph_dim_1", 
                                                                y = "prin_graph_dim_2", label = "branch_point_idx"), 
                                                     size = 4, color = "white", na.rm = TRUE, branch_point_df)
  }
  if (show_cell_names) {
    g <- g + geom_text(aes(label = sample_name), size = cell_name_size)
  }
  if (show_state_number) {
    g <- g + geom_text(aes(label = sample_state), size = state_number_size)
  }
  g <- g + monocle_theme_opts() + xlab(paste("Component", x)) + 
    ylab(paste("Component", y)) + theme(legend.position = "top", 
                                        legend.key.height = grid::unit(0.35, "in")) + theme(legend.key = element_blank()) + 
    theme(panel.background = element_rect(fill = "white"))
  res=list(plot=g,
         edge_df=edge_df,
         data_df=data_df,
         plot_build=ggplot_build(g))
  return(res)
}
#将自定义函数的环境改变为原函数的环境，不然会报错！
environment(my_plot_cell_trajectory)<-environment(plot_cell_trajectory)
#最后可以重新作图啦
a<-my_plot_cell_trajectory(gulou_sc_epi_luminal_Tacstd2_pos_cds)