MUSIC analysis on LUHMES CROP-seq data

Last updated: 2022-09-20

Checks: 7 0

Knit directory: GSFA_analysis/

This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20220524)

The command set.seed(20220524) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: a8e275a

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version a8e275a. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/

Untracked files:
    Untracked:  Rplots.pdf
    Untracked:  analysis/check_Tcells_datasets.Rmd
    Untracked:  analysis/fscLVM_analysis.Rmd
    Untracked:  analysis/spca_LUHMES_data.Rmd
    Untracked:  analysis/test_seurat.Rmd
    Untracked:  code/gsfa_negctrl_job.sbatch
    Untracked:  code/music_LUHMES_Yifan.R
    Untracked:  code/plotting_functions.R
    Untracked:  code/run_fscLVM_LUHMES_data.R
    Untracked:  code/run_gsfa_2groups_negctrl.R
    Untracked:  code/run_gsfa_negctrl.R
    Untracked:  code/run_music_LUHMES.R
    Untracked:  code/run_music_LUHMES_data.sbatch
    Untracked:  code/run_music_LUHMES_data_20topics.R
    Untracked:  code/run_music_LUHMES_data_20topics.sbatch
    Untracked:  code/run_sceptre_Tcells_data.sbatch
    Untracked:  code/run_sceptre_Tcells_stimulated_data.sbatch
    Untracked:  code/run_sceptre_Tcells_test_data.sbatch
    Untracked:  code/run_sceptre_Tcells_unstimulated_data.sbatch
    Untracked:  code/run_sceptre_permuted_data.sbatch
    Untracked:  code/run_spca_LUHMES.R
    Untracked:  code/run_spca_TCells.R
    Untracked:  code/run_twostep_clustering_LUHMES_data.sbatch
    Untracked:  code/run_twostep_clustering_Tcells_data.sbatch
    Untracked:  code/run_unguided_gsfa_LUHMES.R
    Untracked:  code/run_unguided_gsfa_LUHMES.sbatch
    Untracked:  code/run_unguided_gsfa_Tcells.R
    Untracked:  code/run_unguided_gsfa_Tcells.sbatch
    Untracked:  code/sceptre_LUHMES_data.R
    Untracked:  code/sceptre_Tcells_stimulated_data.R
    Untracked:  code/sceptre_Tcells_unstimulated_data.R
    Untracked:  code/sceptre_permutation_analysis.R
    Untracked:  code/sceptre_permute_analysis.R
    Untracked:  code/seurat_sim_fpr_tpr.R
    Untracked:  code/unguided_GFSA_mixture_normal_prior.cpp

Unstaged changes:
    Modified:   analysis/sceptre_TCells_data.Rmd
    Modified:   analysis/twostep_clustering_LUHMES_data.Rmd
    Modified:   code/run_sceptre_LUHMES_data.R
    Modified:   code/run_sceptre_LUHMES_data.sbatch
    Modified:   code/run_sceptre_LUHMES_permuted_data.R
    Modified:   code/run_sceptre_Tcells_permuted_data.R
    Modified:   code/run_sceptre_cropseq_data.sbatch
    Modified:   code/run_twostep_clustering_LUHMES_data.R
    Modified:   code/sceptre_analysis.R

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/music_LUHMES_data.Rmd) and HTML (docs/music_LUHMES_data.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	a8e275a	kevinlkx	2022-09-20	updated dot plots
html	31d48d1	kevinlkx	2022-09-20	Build site.
Rmd	f02521e	kevinlkx	2022-09-20	updated the dot plots and saved as pdf
html	df4b0c6	kevinlkx	2022-08-11	Build site.
Rmd	b979400	kevinlkx	2022-08-11	fix intro formatting
html	a2297e5	kevinlkx	2022-08-11	Build site.
Rmd	a0b124f	kevinlkx	2022-08-11	added dot plots with Bonferroni adjusted p-value
html	43b3a70	kevinlkx	2022-08-01	Build site.
Rmd	706ef41	kevinlkx	2022-08-01	test 4-6 topics
html	e54a27e	kevinlkx	2022-07-29	Build site.
Rmd	66afbb0	kevinlkx	2022-07-29	update dotplot figure legend
html	2763f33	kevinlkx	2022-07-29	Build site.
Rmd	f5bca8a	kevinlkx	2022-07-29	added dotplots with MUSIC empirical TPD FDR
html	223426b	kevinlkx	2022-07-25	Build site.
Rmd	3cbce15	kevinlkx	2022-07-25	update music results for LUHMES data
html	5200b77	kevinlkx	2022-07-19	Build site.
Rmd	358270c	kevinlkx	2022-07-19	run MUSIC on LUHMES data

slurm setting

sinteractive --partition=broadwl --account=pi-xinhe --mem=30G --time=10:00:00 --cpus-per-task=10

MUSIC website: https://github.com/bm2-lab/MUSIC

Scripts for running the analysis:

R script.
sbatch script.

cd /project2/xinhe/kevinluo/GSFA/music_analysis/log

sbatch --mem=50G --cpus-per-task=10 ~/projects/GSFA_analysis/code/run_music_LUHMES_data.sbatch

Running time check

cd /project2/xinhe/kevinluo/GSFA/music_analysis/log

sbatch --mem=50G --cpus-per-task=10 ~/projects/GSFA_analysis/code/run_music_LUHMES_data_20topics.sbatch

About the data sets

CROP-seq datasets: /project2/xinhe/yifan/Factor_analysis/shared_data/. The data are Seurat objects, with raw gene counts stored in obj@assays$RNA@counts, and cell meta data stored in obj@meta.data. Normalized and scaled data used for GSFA are stored in obj@assays$RNA@scale.data, the rownames of which are the 6k genes used for GSFA.

Load packages

dyn.load('/software/geos-3.7.0-el7-x86_64/lib64/libgeos_c.so') # attach the geos lib for Seurat
suppressPackageStartupMessages(library(data.table))
suppressPackageStartupMessages(library(Seurat))
suppressPackageStartupMessages(library(MUSIC))
suppressPackageStartupMessages(library(ComplexHeatmap))
suppressPackageStartupMessages(library(ggplot2))
theme_set(theme_bw() + theme(plot.title = element_text(size = 14, hjust = 0.5),
                             axis.title = element_text(size = 14),
                             axis.text = element_text(size = 13),
                             legend.title = element_text(size = 13),
                             legend.text = element_text(size = 12),
                             panel.grid.minor = element_blank())
)

source("code/plotting_functions.R")

functions

## 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)
}

Set directories

data_dir <- "/project2/xinhe/yifan/Factor_analysis/LUHMES/"
res_dir <- "/project2/xinhe/kevinluo/GSFA/music_analysis/LUHMES"
dir.create(file.path(res_dir,"/music_output"), recursive = TRUE, showWarnings = FALSE)

setwd(res_dir)

Run MUSIC

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)]

cat("Dimension of expression profile matrix: \n")
dim(expression_profile)

targets <- names(combined_obj@meta.data)[4:18]
targets[targets == "Nontargeting"] <- "CTRL"
cat("Targets: \n")
print(targets)

perturb_information <- apply(combined_obj@meta.data[4:18], 1,
                             function(x){ targets[which(x > 0)] })

Dimension of expression profile matrix: 
[1] 33694  8708
Targets: 
 [1] "ADNP"   "ARID1B" "ASH1L"  "CHD2"   "CHD8"   "CTNND2" "DYRK1A" "HDAC5" 
 [9] "MECP2"  "MYT1L"  "CTRL"   "POGZ"   "PTEN"   "RELN"   "SETD5"

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

Obtain highly dispersion differentially expressed genes.

crop_seq_filtered <- readRDS("music_output/music_filtered.merged.rds")
dim(crop_seq_filtered$expression_profile)
length(crop_seq_filtered$perturb_information)

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")

crop_seq_vargene <- readRDS("music_output/music_vargene.merged.rds")
dim(crop_seq_vargene$expression_profile)

crop_seq_vargene$expression_profile[1:5,1:5]
length(crop_seq_vargene$perturb_information)

get topics.

## 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")

## try fewer numbers of topics
system.time(
  topic_5 <- Get_topics(crop_seq_vargene$expression_profile,
                        crop_seq_vargene$perturb_information,
                        topic_number = 4))
saveRDS(topic_5, "music_output/music_merged_4_topics.rds")

system.time(
  topic_6 <- Get_topics(crop_seq_vargene$expression_profile,
                        crop_seq_vargene$perturb_information,
                        topic_number = 6))
saveRDS(topic_6, "music_output/music_merged_6_topics.rds")

3. Pick the number of topics

topic_1 <- readRDS("music_output/music_merged_4_topics.rds")
topic_2 <- readRDS("music_output/music_merged_5_topics.rds")
topic_3 <- readRDS("music_output/music_merged_6_topics.rds")
topic_4 <- readRDS("music_output/music_merged_10_topics.rds")
topic_5 <- readRDS("music_output/music_merged_15_topics.rds")
topic_6 <- readRDS("music_output/music_merged_20_topics.rds")

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]]
topic_model_list$models[[5]] <- topic_5$models[[1]]
topic_model_list$models[[6]] <- topic_6$models[[1]]

optimalModel <- Select_topic_number(topic_model_list$models,
                                    plot = T,
                                    plot_path = "music_output/select_topic_number_4to6to20.pdf")

Summarize the results

Summarize the results under 20 topics to be comparable to GSFA

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(distri_diff, "music_output/distri_diff.rds")

distri_diff <- readRDS("music_output/distri_diff.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()

Overall perturbation effect ranking list.

distri_diff <- readRDS(file.path(res_dir, "music_output/distri_diff.rds"))

rank_overall_result <- Rank_overall(distri_diff)
print(rank_overall_result)
# saveRDS(rank_overall_result, "music_output/rank_overall_result.rds")

   perturbation ranking    Score off_target
1         ASH1L       1 85.91013       none
2         SETD5       2 79.56861       none
3        ARID1B       3 77.59034       none
4          CHD8       4 70.16751       none
5          PTEN       5 66.48328       none
6          POGZ       6 64.12727       none
7         MECP2       7 62.10662       none
8          RELN       8 59.26017       none
9          CHD2       9 55.85265       none
10        MYT1L      10 39.67237       none
11         ADNP      11 39.06289       none
12       DYRK1A      12 38.98911       none
13        HDAC5      13 38.48395       none

topic-specific ranking list.

rank_topic_specific_result <- Rank_specific(distri_diff)
head(rank_topic_specific_result, 10)
# saveRDS(rank_topic_specific_result, "music_output/rank_topic_specific_result.rds")

    topic perturbation ranking
1  Topic1       ARID1B       1
2  Topic1         CHD8       2
3  Topic1        SETD5       3
4  Topic1        MECP2       4
5  Topic1         POGZ       5
6  Topic1         PTEN       6
7  Topic1        HDAC5       7
8  Topic1       DYRK1A       8
9  Topic1         ADNP       9
10 Topic1         RELN      10

Perturbation correlation.

perturb_cor <- Correlation_perturbation(distri_diff,
                                        cutoff = 0.5, gene = "all", plot = T,
                                        plot_path = file.path(res_dir, "music_output/correlation_network_20_topics.pdf"))

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

   Perturbation_1 Perturbation_2 Correlation
2          ARID1B           ADNP -0.36642339
3           ASH1L           ADNP  0.60467107
16          ASH1L         ARID1B  0.12772609
30           CHD2          ASH1L  0.91852307
4            CHD2           ADNP  0.81036468
17           CHD2         ARID1B  0.03762586
18           CHD8         ARID1B  0.89232727
5            CHD8           ADNP -0.26907399
31           CHD8          ASH1L  0.25264119
44           CHD8           CHD2  0.13836576

Adaptation to the code to generate calibrated empirical TPD scores

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(file.path(res_dir, "music_output/music_merged_20_topics_ttest_summary.rds"))
summary_df$topic <- gsub("_", " ", summary_df$topic)
summary_df$topic <- factor(summary_df$topic, levels = paste("Topic", 1:length(unique(summary_df$topic))))

summary_df$fdr <- p.adjust(summary_df$empirical_pval, method = "BH")
summary_df$bonferroni_adj <- p.adjust(summary_df$empirical_pval, method = "bonferroni")

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

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

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

bonferroni_mat <- dcast(summary_df %>% dplyr::select(knockout, topic, bonferroni_adj), knockout ~ topic)
rownames(bonferroni_mat) <- bonferroni_mat$knockout
bonferroni_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)

Version	Author	Date
31d48d1	kevinlkx	2022-09-20
2763f33	kevinlkx	2022-07-29

# dev.off()

# pdf("music_output/music_merged_20_topics_empirical_tstats_fdr_dotplot.pdf",
#     width = 12, height = 8)
KO_names <- rownames(fdr_mat)
dotplot_effectsize(t(effect_mat), t(fdr_mat),
                   reorder_markers = c(KO_names[KO_names!="CTRL"], "CTRL"),
                   color_lgd_title = "MUSIC T statistics",
                   size_lgd_title = "FDR",
                   max_score = 20,
                   min_score = -20,
                   by_score = 10) + 
  coord_flip() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1))

Version	Author	Date
31d48d1	kevinlkx	2022-09-20
43b3a70	kevinlkx	2022-08-01
e54a27e	kevinlkx	2022-07-29
2763f33	kevinlkx	2022-07-29

# dev.off()

KO_names <- rownames(bonferroni_mat)
dotplot_effectsize(t(effect_mat), t(bonferroni_mat),
                   reorder_markers = c(KO_names[KO_names!="CTRL"], "CTRL"),
                   color_lgd_title = "MUSIC T statistics",
                   size_lgd_title = "Bonferroni\nadjusted p-value",
                   max_score = 20,
                   min_score = -20,
                   by_score = 10) + 
  coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1))

Version	Author	Date
31d48d1	kevinlkx	2022-09-20
a2297e5	kevinlkx	2022-08-11

Summarize the results using the optimal number of topics selected by the score

Pick optimal number of topics

topic_1 <- readRDS("music_output/music_merged_4_topics.rds")
topic_2 <- readRDS("music_output/music_merged_5_topics.rds")
topic_3 <- readRDS("music_output/music_merged_6_topics.rds")

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]]

optimalModel <- Select_topic_number(topic_model_list$models,
                                    plot = T,
                                    plot_path = "music_output/select_topic_number_4to6.pdf")

optimalModel
saveRDS(optimalModel, "music_output/optimalModel_4_topics.rds")