WhySingleCell/Figure_1ab.R at main · PayneLab/WhySingleCell · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
library(tidyverse)
library(viridis)
library(ragg)


## Set seed for replication
set.seed(77)

#### PART A ####
## Generate random clusters in a bifurcation structure
points_per_cluster <- 50

### Step 1: generate cluster 1 (it will be the left-most cluster, centered on the y-axis)
x_values <- rnorm(points_per_cluster, mean = 10, sd = 2)
y_values <- rnorm(points_per_cluster, mean = 15, sd = 2)

cluster_1 <- tibble(x_values, y_values) %>%
  mutate(cluster = "Cluster 1")

### Step 2: generate cluster 2 (it will be to the right of cluster 1, centered on the y-axis)
x_values <- rnorm(points_per_cluster, mean = 16, sd = 2)
y_values <- rnorm(points_per_cluster, mean = 15, sd = 2)

cluster_2 <- tibble(x_values, y_values) %>%
  mutate(cluster = "Cluster 2")

### Step 3: generate cluster 3 (it will be to the right of cluster 2, on the top end of the y-axis)
x_values <- rnorm(points_per_cluster, mean = 22, sd = 2)
y_values <- rnorm(points_per_cluster, mean = 20, sd = 2)

cluster_3 <- tibble(x_values, y_values) %>%
  mutate(cluster = "Cluster 3")

### Step 4: generate cluster 4 (it will be to the right of cluster 2, on the bottom end of the y-axis)
x_values <- rnorm(points_per_cluster, mean = 22, sd = 2)
y_values <- rnorm(points_per_cluster, mean = 10, sd = 2)

cluster_4 <- tibble(x_values, y_values) %>%
  mutate(cluster = "Cluster 4")

### Step 5: combine the data points into one table
clusters <- bind_rows(cluster_1, cluster_2, cluster_3, cluster_4)

### Step 6: create data points for a line plot to represent a trajectory
###   overlayed on top of the scatterplot
trajectory_points <- clusters %>%
  group_by(cluster) %>%
  summarize(x_values = mean(x_values), y_values = mean(y_values))

trajectory_line1 <- trajectory_points %>%
  filter(cluster != "Cluster 4")

trajectory_line2 <- trajectory_points %>%
  filter(cluster == "Cluster 2" | cluster == "Cluster 4")

### Step 7: plot the clusters
trajectory_plot <- ggplot() +
  geom_point(clusters, mapping = aes(x = x_values, y = y_values, color = cluster)) +
  scale_color_viridis(discrete = TRUE) +
  # geom_point(clusters, mapping = aes(x = x_values, y = y_values, color = "red")) +
  geom_point(trajectory_points, mapping = aes(x_values, y_values), color = "black", size = 2.5) +
  geom_line(trajectory_line1, mapping = aes(x_values, y_values), color = "black") +
  geom_line(trajectory_line2, mapping = aes(x_values, y_values), color = "black") +
  theme_bw() +
  labs(x = "X", y = "Y", title = "a)") +
  theme(legend.position = "none",
        axis.text.x = element_blank(), axis.ticks.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks.y = element_blank(),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())
  # theme(legend.position = "none",
  #       axis.text.x = element_blank(), axis.ticks.x = element_blank(),
  #       axis.text.y = element_blank(), axis.ticks.y = element_blank())


#### PART B ####
## Split full dataset above into replicates
### Step 1: Randomly assigne a replicate number to each row in the clusters dataframe
replicates <- sample(rep(1:4, each = points_per_cluster))

clusters <- clusters %>%
  mutate(replicate = replicates)

### Step 2: create data points for a line plot to represent a trajectory
###   overlayed on top of the scatterplot
trajectory_points <- clusters %>%
  group_by(replicate, cluster) %>%
  summarize(x_values = mean(x_values), y_values = mean(y_values))

trajectory_line1 <- trajectory_points %>%
  filter(ifelse(replicate != 2,
                cluster != "Cluster 4",
                cluster != "Cluster 3"))

trajectory_line2 <- trajectory_points %>%
  filter(ifelse(replicate != 2,
                cluster == "Cluster 2" | cluster == "Cluster 4",
                cluster == "Cluster 3" | cluster == "Cluster 4"))

### Spet 3: Plot the replicate trajectories
replicate_plot <- ggplot() +
  geom_point(clusters, mapping = aes(x = x_values, y = y_values, color = cluster)) +
  scale_color_viridis(discrete = TRUE) +
  geom_point(trajectory_points, mapping = aes(x_values, y_values), color = "black", size = 2.5) +
  geom_line(trajectory_line1, mapping = aes(x_values, y_values), color = "black") +
  geom_line(trajectory_line2, mapping = aes(x_values, y_values), color = "black") +
  facet_wrap(~replicate) +
  theme_bw() +
  labs(x = "X", y = "Y", title = "b)") +
  theme(legend.position = "none",
        axis.text.x = element_blank(), axis.ticks.x = element_blank(),
        axis.text.y = element_blank(), axis.ticks.y = element_blank(),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())
  # theme(legend.position = "none",
  #       axis.text.x = element_blank(), axis.ticks.x = element_blank(),
  #       axis.text.y = element_blank(), axis.ticks.y = element_blank())


## Save the plots
agg_png("Figure_1a.png", width = 6, height = 6, units = "in", res = 600, scaling = 1.5)
trajectory_plot
dev.off()

agg_png("Figure_1b.png", width = 6, height = 6, units = "in", res = 600, scaling = 1.5)
replicate_plot
dev.off()