simulate_het_for_gRNA <- function(
mtDNA_data,
gRNA_name,
n_perms = 5000,
n_sim_points_to_plot = 2000,
col_pal = NULL
) {
# Filter KO and Control data
KO_df <- mtDNA_data %>%
dplyr::filter(top_guide_final_enrich_combined == gRNA_name) %>%
dplyr::select(Heteroplasmy_all_5024, mtDNA_depth, Depth_all_5024) %>%
dplyr::filter(
Depth_all_5024 >= 0,
Heteroplasmy_all_5024 >= 0,
Heteroplasmy_all_5024 <= 1) %>%
drop_na() %>%
dplyr::mutate(group = gRNA_name)
control_df <- mtDNA_data %>%
dplyr::filter(group == "Control") %>%
dplyr::select(Heteroplasmy_all_5024, mtDNA_depth, Depth_all_5024) %>%
dplyr::filter(
Depth_all_5024 >= 0,
Heteroplasmy_all_5024 >= 0,
Heteroplasmy_all_5024 <= 1) %>%
drop_na() %>%
dplyr::mutate(group = "Control")
# Combine KO and Control data
sim_input_df <- bind_rows(KO_df, control_df)
# Get observed stats
sim_input_stats <- sim_input_df %>%
group_by(group) %>%
dplyr::summarise(
n_cells = n(),
mean_depth = mean(Depth_all_5024, na.rm = TRUE),
var_het = var(Heteroplasmy_all_5024, na.rm = TRUE)
)
obs_ko_var <- sim_input_stats %>% dplyr::filter(group == gRNA_name) %>% dplyr::pull(var_het)
obs_control_var <- sim_input_stats %>% dplyr::filter(group == "Control") %>% dplyr::pull(var_het)
all_sim_var <- numeric(n_perms)
all_sim_hets_list <- vector("list", n_perms)
all_input_hets_list <- vector("list", n_perms)
# Iterate through each permutation
for (perm in seq_len(n_perms)) {
# Sample cells from control group equal to the number of cells in KO group
control_sample <- sim_input_df %>%
dplyr::filter(group == "Control") %>%
sample_n(
min(nrow(control_df), nrow(KO_df)),
replace = FALSE) %>%
dplyr::mutate(group = "Control")
# Sample mtDNA depths from KO group
sampled_depths <- sample(KO_df$Depth_all_5024,
size = nrow(control_sample),
replace = TRUE)
# Sample heteroplasmy levels from control group
sampled_hets <- control_sample$Heteroplasmy_all_5024
all_input_hets_list[[perm]] <- sampled_hets
# Simulate mtDNA copy number
sim_alt_counts <- numeric(nrow(control_sample))
sim_hets <- rep(NA_real_, nrow(control_sample))
valid_indices <- which(sampled_depths > 0)
# Simulate alt counts for valid indices
# where sampled_depths > 0
if (length(valid_indices) > 0) {
sim_alt_counts[valid_indices] <- rbinom(
n = length(valid_indices),
size = sampled_depths[valid_indices],
prob = sampled_hets[valid_indices]
)
# Calculate heteroplasmy for valid indices
sim_hets[valid_indices] <- sim_alt_counts[valid_indices] / sampled_depths[valid_indices]
}
all_sim_var[perm] <- var(sim_hets, na.rm = TRUE)
all_sim_hets_list[[perm]] <- sim_hets[!is.na(sim_hets)]
}
# Flatten the list of simulated heteroplasmies
all_sim_hets_flat <- na.omit(unlist(all_sim_hets_list))
# Observed KO stats
obs_ko_homoplasmy_prop <- mean(KO_df$Heteroplasmy_all_5024 %in% c(0,1), na.rm = TRUE)
obs_ko_mean_het <- mean(KO_df$Heteroplasmy_all_5024, na.rm = TRUE)
# Simulated stats
sim_ko_homoplasmy_prop <- mean(all_sim_hets_flat %in% c(0,1), na.rm = TRUE)
sim_ko_mean_het <- mean(all_sim_hets_flat, na.rm = TRUE)
sim_ko_var <- var(all_sim_hets_flat, na.rm = TRUE)
# Empirical p-values
# Two-tailed: proportion of simulated variances as extreme or more extreme than observed
mean_sim_var <- mean(all_sim_var)
p_val_two_tailed <- mean(abs(all_sim_var - mean_sim_var) >= abs(obs_ko_var - mean_sim_var))
plot_df_hets <- data.frame(
Heteroplasmy = c(
control_df$Heteroplasmy_all_5024,
all_sim_hets_flat,
KO_df$Heteroplasmy_all_5024
),
Group = factor(c(
rep("Control", nrow(control_df)),
rep("Simulated", length(all_sim_hets_flat)),
rep(gRNA_name, nrow(KO_df))
), levels = c("Control", "Simulated", gRNA_name))
)
# Subsample points to avoid cluttering
# Set the number of points to plot for the jitter layer
n_sim_points <- min(length(all_sim_hets_flat), n_sim_points_to_plot)
# Filter the main dataframe, subsample only the 'Simulated' group
jitter_data <- plot_df_hets %>%
group_by(Group) %>%
# If the group is 'Simulated', sample n_sim_points, otherwise take all rows
dplyr::filter(
if (cur_group()$Group == "Simulated") row_number() %in% sample(row_number(), n_sim_points) else TRUE
) %>%
ungroup()
# Set colors for Control, Simulated, and gRNA
if (!is.null(col_pal)) {
group_colors <- c(
"Control" = scales::alpha(shades::saturation(col_pal["NT"], 0.8), 0.6),
"Simulated" = scales::alpha(shades::saturation(col_pal[gRNA_name], 0.4), 0.6),
gRNA_name = scales::alpha(shades::saturation(col_pal[gRNA_name], 0.8), 0.6)
)
names(group_colors)[3] <- gRNA_name
} else {
group_colors <- c(
"Control" = "#2CA02CFF",
"Simulated" = "#1F77B4FF",
gRNA_name = "#D62728FF"
)
names(group_colors)[3] <- gRNA_name
}
het_dist_sim <- ggplot(
plot_df_hets,
aes(x = Group, y = Heteroplasmy, fill = Group)
) +
geom_violin(
data = jitter_data,
width = 1.2,
alpha = 0.3,
color = NA
) +
geom_jitter(
data = jitter_data,
width = 0.15,
size = 1.2,
alpha = 0.3,
aes(color = Group)
) +
geom_boxplot(
data = jitter_data,
width = 0.2,
outlier.shape = NA,
alpha = 0.6,
position = "identity"
) +
scale_fill_manual(values = group_colors) +
scale_color_manual(values = group_colors) +
theme_minimal(base_size = 14) +
labs(
title = paste("Distribution of Heteroplasmies (", gRNA_name, " KO Simulation)", sep = ""),
x = NULL,
y = "Heteroplasmy (m.5024)"
) +
theme(legend.position = "none")
return(list(
plot = het_dist_sim,
sim_input_stats = sim_input_stats,
all_sim_var = all_sim_var,
all_sim_hets_list = all_sim_hets_list,
obs_ko_mean_het = obs_ko_mean_het,
obs_ko_var = obs_ko_var,
obs_ko_homoplasmy_prop = obs_ko_homoplasmy_prop,
sim_ko_mean_het = sim_ko_mean_het,
sim_ko_var = sim_ko_var,
sim_ko_homoplasmy_prop = sim_ko_homoplasmy_prop,
p_value_two_tailed = p_val_two_tailed
))
}
