Rename, build man, build covid vignette.

R/pppr_output_plot.R

@@ -0,0 +1,117 @@
+#' PPR OUTPUT PLOT
+#'
+#' @description
+#' Plots the predicted positon of your sample.
+#'
+#' @param sample.gss A GSS_OBJECT of the sample you wish to plot
+#' @param col The colour that you'd like
+#' @param overlay set to TRUE if you would like this plot to overlay a previous plot.
+#' @param label string that you would like to assign as the label to the line.
+#' @param genes_of_interest The names of any genes that you'd like to include the switching point of.
+#' @param switching_genes a matrix containing the switching gene information as produced by GeneSwitches.
+#'
+#' @return nice plot highlighting the probable position of your sample on your trajectory.
+#' @export
+#'
+ppr_output_plot <- function(sample.gss, col = "red", overlay = FALSE, label = "sample name", genes_of_interest = NULL, switching_genes = reference.sg){
+
+ if (!overlay) {
+  plot(x = 1:100,
+       y = (sample.gss$sample_flat/max(sample.gss$sample_flat)*100),
+       ylim = c(0,110),
+       xlim = c(0,100),
+       pch = 20,
+       cex = 0.8,
+       col = col,
+       type = "l",
+       xlab = "Pseudo-Time Index",
+       ylab = "GSS Score",
+       main = paste("Predicted Positions"))
+
+  segments(which.max(colSums(sample.gss$sample_flat)),
+           -3.9,
+           which.max(colSums(sample.gss$sample_flat)),
+           (sample.gss$sample_flat[which.max(colSums(sample.gss$sample_flat))]/max(sample.gss$sample_flat)*100),
+           lwd = 1,
+           lty = 2,
+           col = col)
+
+  text(x = which.max(colSums(sample.gss$sample_flat)),
+       y = (sample.gss$sample_flat[which.max(colSums(sample.gss$sample_flat))]/max(sample.gss$sample_flat)*100),
+       labels = label,
+       col = col,
+       pos = 3)
+
+  if (length(genes_of_interest) > 0) {
+    for (gene_name in genes_of_interest) {
+      segments(switching_genes[gene_name,"switch_at_timeidx"],
+               -3.9,
+               switching_genes[gene_name,"switch_at_timeidx"],
+               -0.5,
+               lwd = 1,
+               lty = 2)
+
+      text(x = switching_genes[gene_name,"switch_at_timeidx"] + 3,
+           y = 1,
+           labels = gene_name,
+           srt = -20,
+           cex = 0.86,
+           pos = 2)
+
+    }
+  }
+
+
+
+  } else {
+    lines(x = 1:100,
+         y = (sample.gss$sample_flat/max(sample.gss$sample_flat)*100),
+         ylim = c(0,110),
+         xlim = c(0,100),
+         pch = 20,
+         cex = 0.8,
+         col = col,
+         type = "l")
+
+    segments(which.max(colSums(sample.gss$sample_flat)),
+             -3.9,
+             which.max(colSums(sample.gss$sample_flat)),
+             (sample.gss$sample_flat[which.max(colSums(sample.gss$sample_flat))]/max(sample.gss$sample_flat)*100),
+             lwd = 1,
+             lty = 2,
+             col = col)
+
+    text(x = which.max(colSums(sample.gss$sample_flat)),
+         y = (sample.gss$sample_flat[which.max(colSums(sample.gss$sample_flat))]/max(sample.gss$sample_flat)*100),
+         labels = label,
+         col = col,
+         pos = 3)
+
+    if (length(genes_of_interest) > 0) {
+      for (gene_name in genes_of_interest) {
+        segments(switching_genes[gene_name,"switch_at_timeidx"],
+                 -3.9,
+                 switching_genes[gene_name,"switch_at_timeidx"],
+                 -0.5,
+                 lwd = 1,
+                 lty = 2)
+
+        text(x = switching_genes[gene_name,"switch_at_timeidx"] + 3,
+             y = 1,
+             labels = gene_name,
+             srt = -20,
+             cex = 0.86,
+             pos = 2)
+
+   }
+  }
+ }
+}
+
+
+
+
+
+
+
+

R/pppr_predict_position.R

@@ -0,0 +1,88 @@
+#' @title PPPR Predict Position
+#'
+#' @description
+#' Produces an estimate for the position on trajectory of each gene in each cell of a sample.
+#'     This can be later aggregated to estimate the position of the sample along the trajectory.
+#'
+#' @param reduced_binary_counts_matrix a matrix of your samples binary gene expression.
+#' @param reference.sg Genes which switch through the trajectory as identified by GeneSwitches.
+#'
+#' @return A list of matrices: A matrix for each cell where the columns represent progress through a trajectory,
+#'     and the rows represent genes, values indicate a likely position of the cell upon the trajectory based that genes bianrized expression.
+#' @export
+#'
+ppr_predict_position <- function(reduced_binary_counts_matrix,reference.sg) {
+
+  ## The final output will be a ppr_obj (list) comprised of 3 objects.
+  # Make the list of length 3, and name the objects
+  ppr_obj <- vector("list", 3)
+  names(ppr_obj) <- c("genomic_expression_traces","cells_flat","sample_flat")
+  # Assign the ppr_obj class attribute to the list
+  class(ppr_obj) <- "PPR_OBJECT"
+
+  ## Reorder reference.sg
+  # as the code relies on the rownames and idicies of the genes in reduced_binary_counts_matrix and reference.sg matching.
+  reference.sg <- reference.sg[rownames(reduced_binary_counts_matrix),]
+
+  ## Input values:
+  number_of_cells <- dim(reduced_binary_counts_matrix)[2]
+  # Should we get the number of switching genes from here? or dim(reduced_binary_counts_matrix)[1]
+  # It is possible for dim(reduced_binary_counts_matrix)[1] <  dim(reference.sg)[1]
+  number_of_switching_genes <- dim(reference.sg)[1]
+
+
+  # Different datasets store this information in different columns.. watch out.
+  # Maybe this should be an input to the function?
+  # as.numeric is needed/notneeded depending on dataset.
+  switching_time <- as.numeric(reference.sg$switch_at_timeidx)
+  switching_direction <- reference.sg$direction
+
+  # Building the genomic_expression_traces list. (faster than building it dynamically.)
+  all_patients_cells_scored <- vector("list", number_of_cells)
+  names(all_patients_cells_scored) <- colnames(reduced_binary_counts_matrix)
+
+  # Loop through all cells making matrices for each,
+  #which represent likely position of a cell on a trajectory based on the expression of each gene.
+  for (c in 1:number_of_cells) {
+    # Build the matrix of 0's which has genes as rows and pseudotime indecies as columns.
+   genomic_expression_mat <- matrix(0, nrow = number_of_switching_genes, ncol = 100)
+   rownames(genomic_expression_mat) <- rownames(reduced_binary_counts_matrix)
+    binarized_gene_expression_for_cell_c <- reduced_binary_counts_matrix[, c]
+
+    up_indices <- which(binarized_gene_expression_for_cell_c == 1 & switching_direction == "up")
+    down_indices <- which(binarized_gene_expression_for_cell_c == 0 & switching_direction == "down")
+    not_up_indices <- which(binarized_gene_expression_for_cell_c == 0 & switching_direction == "up")
+    not_down_indices <- which(binarized_gene_expression_for_cell_c == 1 & switching_direction == "down")
+
+    for (i in up_indices) {
+      genomic_expression_mat[i, switching_time[i]:100] <- 1
+    }
+
+    for (i in down_indices) {
+      genomic_expression_mat[i, switching_time[i]:100] <- 1
+    }
+
+    for (i in not_up_indices) {
+      genomic_expression_mat[i, 1:switching_time[i]] <- 1
+    }
+
+    for (i in not_down_indices) {
+      genomic_expression_mat[i, 1:switching_time[i]] <- 1
+    }
+
+    all_patients_cells_scored[[c]] <- genomic_expression_mat
+  }
+
+  ppr_obj$genomic_expression_traces <- all_patients_cells_scored
+
+### GENOMIC EXPRESSION TRACES CREATED
+# Now flatten:
+
+# Use lapply to calculate column sums for each matrix
+ppr_obj$cells_flat <- do.call(rbind, lapply(all_patients_cells_scored, colSums))
+rownames(ppr_obj$cells_flat) <- names(all_patients_cells_scored)
+# Combine each cells column sums into a single flat matrix.
+ppr_obj$sample_flat <- matrix(colSums(ppr_obj$cells_flat), nrow = 1)
+
+return(ppr_obj)
+}

R/pppr_racinglines_plot.R

@@ -0,0 +1,113 @@
+#' @title Identify and Visualise each cells position.
+#'
+#' @description
+#' Produces a plot for each cell which helps visualize how PathPinPointR is predicting the cells position.
+#'
+#' @param reference.sg A selection of switching genes which are evenly distributed through pseudo-time.
+#' @param lines logical, do you want to include lines to inidicate the predicted position of a cell
+#' @param reference_reduced a matrix of your sample's binary gene expression.
+#' @param cell_id The index or name of the cell of interest
+#'
+#' @return Timeline plot of selected cell
+#' @import ggplot2
+#' @import ggrepel
+#'
+#'
+#' @export
+#'
+
+ppr_timeline_plot <- function(reference.sg, lines = TRUE, reference_reduced = NULL, cell_id = 1) {
+
+  # Convert reference.sg to a data frame
+  reference.sg <- as.data.frame(reference.sg)
+
+  # Add a new column direction_num and set it to -1
+  reference.sg$direction_num <- -1
+
+  # If "up" is present in the direction column, set direction_num to 1
+  if ("up" %in% reference.sg$direction) {
+    reference.sg[reference.sg$direction == "up", ]$direction_num <- 1
+  }
+
+  # Generate pseudotime data based on the specified parameters
+    timeidx_df <- data.frame(timeidx_range = c(0, 25, 50, 75, 100), timeidx_labs = c(0, 25, 50, 75, 100))
+
+
+  # Create the initial ggplot object with x and y aesthetics, color, and labels
+  timeline_plot <- ggplot(reference.sg, aes(x = switch_at_timeidx, y = pseudoR2s * direction_num, label = rownames(reference.sg))) +
+    geom_point(size = 1) + xlab("Time-Index") + ylab("Quality of fitting (R^2)")
+
+  # Add the classic theme to the plot
+  timeline_plot <- timeline_plot + theme_classic()
+
+  # Add a horizontal black line for the timeline
+  timeline_plot <- timeline_plot + geom_hline(yintercept = 0, color = "black", linewidth = 0.6)
+
+  # Add labels for pseudotime on the plot
+  timeline_plot <- timeline_plot + geom_label(data = timeidx_df, aes(x = timeidx_range, y = 0, label = timeidx_labs),
+                                    size = (3), color = "black")
+
+  # Add text labels with repulsion to avoid overlap
+  timeline_plot <- timeline_plot + geom_text_repel(aes(x = switch_at_timeidx, y = pseudoR2s * direction_num, label = rownames(reference.sg)),
+                                         size = 3, show.legend = FALSE)
+
+  # Customize the theme and legend appearance
+  timeline_plot <- timeline_plot + theme(legend.position = "bottom", legend.title = element_blank(), legend.key.size = unit(10, "pt"),
+                               text = element_text(size = 12))
+
+
+if (lines) {
+
+  if (is.character(cell_id)) {
+    cell_idx <- which(colnames(reference_reduced) == cell_id)
+  } else if (is.numeric(cell_id)) {
+    # Convert numeric to integer if needed
+    cell_idx <- as.integer(cell_id)
+  } else {
+    stop("cell_id must be character or numeric.")
+  }
+
+  ## Reorder reference.sg
+  # as the code relies on the rownames and idicies of the genes in reference_reduced and reference.sg matching.
+  reference_reduced <- reference_reduced[rownames(reference.sg),]
+
+  # Check if row names in reference.sg match those in reference_reduced
+  if (!setequal(rownames(reference.sg), rownames(reference_reduced))) {
+    stop("Row names in reference.sg do not match those in reference_reduced")
+  }
+
+  #loop through all of the genes in reference.sg.
+    for (g in 1:dim(reference.sg)[1]) {
+      # if G is NOT expressed in  C, and the switch is UP  , then draw the line to the right.
+      if ((reference_reduced[rownames(reference.sg)[g], cell_idx] == 0) && (reference.sg$direction[g] == "up")) {
+        timeline_plot <- timeline_plot + geom_segment(aes_string(x = 0, xend = reference.sg$switch_at_timeidx[g], y = reference.sg$pseudoR2s[g], yend = reference.sg$pseudoR2s[g]),
+                                                      color = "blue", linewidth = 0.6)
+      }
+      # if G is expressed in      c, and the switch is UP  , then draw the line to the right.
+      if ((reference_reduced[rownames(reference.sg)[g], cell_idx] == 1) && (reference.sg$direction[g] == "up")) {
+        timeline_plot <- timeline_plot + geom_segment(aes_string(x = reference.sg$switch_at_timeidx[g], xend = 100, y = reference.sg$pseudoR2s[g], yend = reference.sg$pseudoR2s[g]),
+                                                      color = "blue", linewidth = 0.6)
+      }
+      # if G is NOT expressed in  C, and the switch is Down, then draw the line to the Left.
+      if ((reference_reduced[rownames(reference.sg)[g], cell_idx] == 0) && (reference.sg$direction[g] == "down")) {
+        timeline_plot <- timeline_plot + geom_segment(aes_string(x = reference.sg$switch_at_timeidx[g], xend = 100, y = -reference.sg$pseudoR2s[g], yend = -reference.sg$pseudoR2s[g]),
+                                                      color = "blue", linewidth = 0.6)
+      }
+      # if G is expressed in      C, and the switch is Down, then draw the line to the Left.
+      if ((reference_reduced[rownames(reference.sg)[g], cell_idx] == 1) && (reference.sg$direction[g] == "down")) {
+        timeline_plot <- timeline_plot + geom_segment(aes_string(x = 0, xend = reference.sg$switch_at_timeidx[g], y = -reference.sg$pseudoR2s[g], yend = -reference.sg$pseudoR2s[g]),
+                                                      color = "blue", linewidth = 0.6)
+      }
+    }
+
+# Title the plot with the name of the cell
+timeline_plot <- timeline_plot + ggtitle(colnames(reference_reduced)[cell_idx])
+
+
+}
+
+# Return the final plot
+return(timeline_plot)
+
+}
+