MUSIC analysis on LUHMES CROP-seq data

0. Load input data

feature.names <- data.frame(fread(paste0(data_dir, "GSE142078_raw/GSM4219575_Run1_genes.tsv.gz"),
                                  header = FALSE), stringsAsFactors = FALSE)

# combined_obj <- readRDS("processed_data/seurat_obj.merged_scaled_detect_01.corrected_new.rds")
combined_obj <- readRDS("/project2/xinhe/yifan/Factor_analysis/shared_data/LUHMES_cropseq_data_seurat.rds")
expression_profile <- combined_obj@assays$RNA@counts
rownames(expression_profile) <- feature.names$V2[match(rownames(expression_profile),
                                                       feature.names$V1)]
targets <- names(combined_obj@meta.data)[4:18]
targets[11] <- "CTRL"
perturb_information <- apply(combined_obj@meta.data[4:18], 1,
                             function(x){ targets[which(x > 0)] })

1. Data preprocessing

crop_seq_list <- Input_preprocess(expression_profile, perturb_information)

crop_seq_qc <- Cell_qc(crop_seq_list$expression_profile,
                       crop_seq_list$perturb_information,
                       species = "Hs", plot = F)

crop_seq_imputation <- Data_imputation(crop_seq_qc$expression_profile,
                                       crop_seq_qc$perturb_information,
                                       cpu_num = 10)
saveRDS(crop_seq_imputation, "music_output/music_imputation.merged.rds")

crop_seq_filtered <- Cell_filtering(crop_seq_imputation$expression_profile,
                                    crop_seq_imputation$perturb_information,
                                    cpu_num = 10)
saveRDS(crop_seq_filtered, "music_output/music_filtered.merged.rds")

2. Model building

crop_seq_vargene <- Get_high_varGenes(crop_seq_filtered$expression_profile,
                                      crop_seq_filtered$perturb_information, plot = T)
saveRDS(crop_seq_vargene, "music_output/music_vargene.merged.rds")

## Get_topics() can take up to a few hours to finish, 
## depending on the size of data
system.time(
  topic_1 <- Get_topics(crop_seq_vargene$expression_profile,
                        crop_seq_vargene$perturb_information,
                        topic_number = 5))
saveRDS(topic_1, "music_output/music_merged_5_topics.rds")

system.time(
  topic_2 <- Get_topics(crop_seq_vargene$expression_profile,
                        crop_seq_vargene$perturb_information,
                        topic_number = 10))
saveRDS(topic_2, "music_output/music_merged_10_topics.rds")

system.time(
  topic_3 <- Get_topics(crop_seq_vargene$expression_profile,
                        crop_seq_vargene$perturb_information,
                        topic_number = 15))
saveRDS(topic_3, "music_output/music_merged_15_topics.rds")

system.time(
  topic_4 <- Get_topics(crop_seq_vargene$expression_profile,
                        crop_seq_vargene$perturb_information,
                        topic_number = 20))
saveRDS(topic_4, "music_output/music_merged_20_topics.rds")

3. Pick optimal number of topics

topic_model_list <- list()
topic_model_list$models <- list()
topic_model_list$perturb_information <- topic_1$perturb_information
topic_model_list$models[[1]] <- topic_1$models[[1]]
topic_model_list$models[[2]] <- topic_2$models[[1]]
topic_model_list$models[[3]] <- topic_3$models[[1]]
topic_model_list$models[[4]] <- topic_4$models[[1]]

saveRDS(topic_model_list, "music_output/topic_model_list.rds")

optimalModel <- Select_topic_number(topic_model_list$models,
                                    plot = T,
                                    plot_path = "music_output/select_topic_number_5_to_20.pdf")
# The paper said "the larger the score, the better the selected topic number".
# But we probably are going to set the number to 20 just to be comparable to GSFA.

4. Visualize the results under 20 topics

Gene ontology annotations for top topics

topic_res <- readRDS("music_output/music_merged_20_topics.rds")

topic_func <- Topic_func_anno(topic_res$models[[1]], species = "Hs")
saveRDS(topic_func, "music_output/topic_func.rds")

topic_func <- readRDS("music_output/topic_func.rds")
# pdf("music_output/music_merged_20_topics_GO_annotations.pdf",
#     width = 14, height = 12)
ggplot(topic_func$topic_annotation_result) +
  geom_point(aes(x = Cluster, y = Description,
                 size = Count, color = -log10(qvalue))) +
  scale_color_gradientn(colors = c("blue", "red")) +
  theme_bw() +
  theme(axis.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1))
# dev.off()

Perturbation effect prioritizing

# calculate topic distribution for each cell.
distri_diff <- Diff_topic_distri(topic_res$models[[1]],
                                 topic_res$perturb_information,
                                 plot = T)
saveRDS(topic_func, "music_output/topic_func.rds")

t_D_diff_matrix <- dcast(distri_diff %>% dplyr::select(knockout, variable, t_D_diff),
                         knockout ~ variable)
rownames(t_D_diff_matrix) <- t_D_diff_matrix$knockout
t_D_diff_matrix$knockout <- NULL

# pdf("music_output/music_merged_20_topics_TPD_heatmap.pdf", width = 12, height = 8)
Heatmap(t_D_diff_matrix,
        name = "Topic probability difference (vs ctrl)",
        cluster_rows = T, cluster_columns = T,
        column_names_rot = 45,
        heatmap_legend_param = list(title_gp = gpar(fontsize = 12, fontface = "bold")))
# dev.off()

# calculate the overall perturbation effect ranking list without "offTarget_Info".
rank_overall_result <- Rank_overall(distri_diff)
head(rank_overall_result)
saveRDS(rank_overall_result, "music_output/rank_overall_result.rds")

# calculate the topic-specific ranking list.
rank_topic_specific_result <- Rank_specific(distri_diff)
head(rank_topic_specific_result)
saveRDS(rank_topic_specific_result, "music_output/rank_topic_specific_result.rds")

# calculate the perturbation correlation.
perturb_cor <- Correlation_perturbation(distri_diff,
                                        cutoff = 0.5, gene = "all", plot = T,
                                        plot_path = "music_output/correlation_network_20_topics.pdf")

head(perturb_cor)
saveRDS(perturb_cor, "music_output/perturb_cor.rds")

5. Adaptation to the code to generate calibrated empirical TPD scores

## Adapted over MUSIC's Diff_topic_distri() function
Empirical_topic_prob_diff <- function(model, perturb_information,
                                      permNum = 10^4, seed = 1000){
  require(reshape2)
  require(dplyr)
  require(ComplexHeatmap)
  options(warn = -1)
  prob_mat <- model@gamma
  row.names(prob_mat) <- model@documents
  topicNum <- ncol(prob_mat)
  topicName <- paste0('Topic_', 1:topicNum)
  colnames(prob_mat) <- topicName
  ko_name <- unique(perturb_information)
  prob_df <- data.frame(prob_mat, 
                        samples = rownames(prob_mat),
                        knockout = perturb_information)
  
  prob_df <- melt(prob_df, id = c('samples', 'knockout'), variable.name = "topic")
  
  summary_df <- prob_df %>%
    group_by(knockout, topic) %>%
    summarise(number = sum(value)) %>%
    ungroup() %>%
    group_by(knockout) %>%
    mutate(cellNum = sum(number)) %>%
    ungroup() %>%
    mutate(ratio = number/cellNum)
  
  summary_df$ctrlNum <- rep(summary_df$cellNum[summary_df$knockout == "CTRL"],
                            length(ko_name))
  summary_df$ctrl_ratio <- rep(summary_df$ratio[summary_df$knockout == "CTRL"],
                               length(ko_name))
  summary_df <- summary_df %>% mutate(diff_index = ratio - ctrl_ratio)
  
  test_df <- data.frame(matrix(nrow = length(ko_name) * topicNum, ncol = 5))
  colnames(test_df) <- c("knockout", "topic", "obs_t_stats", "obs_pval", "empirical_pval")
  k <- 1
  for(i in topicName){
    prob_df.topic <- prob_df[prob_df$topic == i, ]
    ctrl_topic <- prob_df.topic$value[prob_df.topic$knockout == "CTRL"]
    ctrl_topic_z <- (ctrl_topic - mean(ctrl_topic)) / sqrt(var(ctrl_topic))
    for(j in ko_name){
      ko_topic <- prob_df.topic$value[prob_df.topic$knockout == j]
      ko_topic_z <- (ko_topic - mean(ctrl_topic)) / sqrt(var(ctrl_topic))
      test_df$knockout[k] <- j
      test_df$topic[k] <- i
      test <- t.test(ko_topic_z, ctrl_topic_z)
      test_df$obs_t_stats[k] <- test$statistic
      test_df$obs_pval[k] <- test$p.value
      k <- k + 1
    }
  }
  
  ## Permutation on the perturbation conditions:
  # permNum <- 10^4
  print(paste0("Performing permutation for ", permNum, " rounds."))
  perm_t_stats <- matrix(0, nrow = nrow(test_df), ncol = permNum)
  set.seed(seed)
  for (perm in 1:permNum){
    perm_prob_df <- data.frame(prob_mat, 
                               samples = rownames(prob_mat),
                               knockout = perturb_information[sample(length(perturb_information))])
    perm_prob_df <- melt(perm_prob_df, id = c('samples', 'knockout'), variable.name = "topic")
    k <- 1
    for(i in topicName){
      perm_prob_df.topic <- perm_prob_df[perm_prob_df$topic == i, ]
      ctrl_topic <- perm_prob_df.topic$value[perm_prob_df.topic$knockout == "CTRL"]
      ctrl_topic_z <- (ctrl_topic - mean(ctrl_topic)) / sqrt(var(ctrl_topic))
      for(j in ko_name){
        ko_topic <- perm_prob_df.topic$value[perm_prob_df.topic$knockout == j]
        ko_topic_z <- (ko_topic - mean(ctrl_topic)) / sqrt(var(ctrl_topic))
        test <- t.test(ko_topic_z, ctrl_topic_z)
        perm_t_stats[k, perm] <- test$statistic
        k <- k + 1
      }
    }
    if (perm %% 1000 == 0){
      print(paste0(perm, " rounds finished."))
    }
  }
  ## Compute two-sided empirical p value:
  for (k in 1:nrow(test_df)){
    test_df$empirical_pval[k] <-
      2 * min(mean(perm_t_stats[k, ] <= test_df$obs_t_stats[k]),
              mean(perm_t_stats[k, ] >= test_df$obs_t_stats[k]))
  }
  test_df <- test_df %>%
    mutate(empirical_pval = ifelse(empirical_pval == 0, 1/permNum, empirical_pval)) %>%
    mutate(empirical_pval = ifelse(empirical_pval > 1, 1, empirical_pval))
  
  summary_df <- inner_join(summary_df, test_df, by = c("knockout", "topic"))
  summary_df <- summary_df %>%
    mutate(polar_log10_pval = ifelse(obs_t_stats > 0, -log10(empirical_pval), log10(empirical_pval)))
  return(summary_df)
}

summary_df <- Empirical_topic_prob_diff(topic_res$models[[1]],
                                        topic_res$perturb_information)
saveRDS(summary_df, "music_output/music_merged_20_topics_ttest_summary.rds")

summary_df <- readRDS("music_output/music_merged_20_topics_ttest_summary.rds")

log10_pval_mat <- dcast(summary_df %>% dplyr::select(knockout, topic, polar_log10_pval),
                        knockout ~ topic)
rownames(log10_pval_mat) <- log10_pval_mat$knockout
log10_pval_mat$knockout <- NULL

pdf("music_output/music_merged_20_topics_empirical_tstats_heatmap.pdf",
    width = 12, height = 8)
ht <- Heatmap(log10_pval_mat,
              name = "Polarized empirical t-test -log10(p-value)\n(KO vs ctrl cell topic probs)",
              col = circlize::colorRamp2(breaks = c(-4, 0, 4), colors = c("blue", "grey90", "red")),
              cluster_rows = T, cluster_columns = T,
              column_names_rot = 45,
              heatmap_legend_param = list(title_gp = gpar(fontsize = 12,
                                                          fontface = "bold")))
draw(ht)
dev.off()

Session information

sessionInfo()