MUSIC analysis on CD8+ T Cell (stimulated and unstimulated) CROP-seq data

Last updated: 2022-07-25

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: 3178265

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 3178265. 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:    .Rhistory
    Ignored:    .Rproj.user/

Untracked files:
    Untracked:  analysis/interpret_gsfa_LUHMES.Rmd
    Untracked:  analysis/spca_LUHMES_data.Rmd
    Untracked:  code/music_LUHMES_Yifan.R
    Untracked:  code/plotting_functions.R
    Untracked:  code/run_music_LUHMES.R
    Untracked:  code/run_music_LUHMES_data.sbatch
    Untracked:  code/run_music_Tcells_stimulated_data.R
    Untracked:  code/run_music_Tcells_stimulated_data.sbatch
    Untracked:  code/run_music_Tcells_unstimulated_data.R
    Untracked:  code/run_music_Tcells_unstimulated_data.sbatch
    Untracked:  code/run_sceptre_LUHMES_data.sbatch
    Untracked:  code/run_sceptre_Tcells_stimulated_data.sbatch
    Untracked:  code/run_sceptre_Tcells_unstimulated_data.sbatch
    Untracked:  code/run_spca_LUHMES.R
    Untracked:  code/run_spca_TCells.R
    Untracked:  code/sceptre_LUHMES_data.R
    Untracked:  code/sceptre_Tcells_stimulated_data.R
    Untracked:  code/sceptre_Tcells_unstimulated_data.R

Unstaged changes:
    Modified:   analysis/music_LUHMES_data.Rmd
    Modified:   analysis/sceptre_LUHMES_data.Rmd
    Modified:   code/run_sceptre_cropseq_data.sbatch
    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_TCells_data.Rmd) and HTML (docs/music_TCells_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	3178265	kevinlkx	2022-07-25	update music results for T cells

slurm setting

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

Scripts for running the analysis:

R script.
sbatch script.

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

sbatch --cpus-per-task=6 --mem=50G ~/projects/GSFA_analysis/code/run_music_Tcells_stimulated_data.sbatch
sbatch --cpus-per-task=6 --mem=40G ~/projects/GSFA_analysis/code/run_music_Tcells_stimulated_data.sbatch

sbatch --partition=xinhe --account=pi-xinhe --cpus-per-task=6 --mem=40G ~/projects/GSFA_analysis/code/run_music_Tcells_stimulated_data.sbatch

sbatch --partition=mstephens --account=pi-mstephens --cpus-per-task=6 --mem=100G ~/projects/GSFA_analysis/code/run_music_Tcells_stimulated_data.sbatch

sbatch --partition=bigmem2 --cpus-per-task=6 --mem=100G ~/projects/GSFA_analysis/code/run_music_Tcells_stimulated_data.sbatch

sbatch --mem=50G --cpus-per-task=5 ~/projects/GSFA_analysis/code/run_music_Tcells_unstimulated_data.sbatch

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

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

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

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

Set directories

data_dir <- "/project2/xinhe/yifan/Factor_analysis/Stimulated_T_Cells/"
dir.create("/project2/xinhe/kevinluo/GSFA/music_analysis/Stimulated_T_Cells", recursive = TRUE, showWarnings = FALSE)

Load the T Cells CROP-seq data

combined_obj <- readRDS('/project2/xinhe/yifan/Factor_analysis/shared_data/TCells_cropseq_data_seurat.rds')

Separate stimulated and unstimulated cells into two data sets, and run those separately.

metadata <- combined_obj@meta.data
metadata[1:5, ]

table(metadata$orig.ident)

stimulated_cells <- rownames(metadata)[which(endsWith(metadata$orig.ident, "S"))]
cat(length(stimulated_cells), "stimulated cells. \n")
unstimulated_cells <- rownames(metadata)[which(endsWith(metadata$orig.ident, "N"))]
cat(length(unstimulated_cells), "unstimulated cells. \n")

                     orig.ident nCount_RNA nFeature_RNA ARID1A BTLA C10orf54
D1S_AAACCTGAGCGATTCT TCells_D1S       8722         2642      0    0        0
D1S_AAACCTGCAAGCCGTC TCells_D1S       7296         2356      0    1        0
D1S_AAACCTGCAGGGCATA TCells_D1S       6467         2154      0    0        0
D1S_AAACCTGGTAAGGGCT TCells_D1S      12151         3121      0    0        0
D1S_AAACCTGGTAGCGCTC TCells_D1S      11170         2964      0    0        1
                     CBLB CD3D CD5 CDKN1B DGKA DGKZ HAVCR2 LAG3 LCP2 MEF2D
D1S_AAACCTGAGCGATTCT    0    0   0      0    0    0      0    0    0     0
D1S_AAACCTGCAAGCCGTC    0    0   0      0    0    0      0    0    0     0
D1S_AAACCTGCAGGGCATA    0    0   0      0    0    0      0    0    0     0
D1S_AAACCTGGTAAGGGCT    1    0   0      0    0    0      0    0    0     0
D1S_AAACCTGGTAGCGCTC    0    0   0      0    0    0      0    0    0     0
                     NonTarget PDCD1 RASA2 SOCS1 STAT6 TCEB2 TMEM222 TNFRSF9
D1S_AAACCTGAGCGATTCT         0     0     0     0     0     1       0       0
D1S_AAACCTGCAAGCCGTC         0     0     0     0     0     0       0       0
D1S_AAACCTGCAGGGCATA         0     0     0     0     1     0       0       0
D1S_AAACCTGGTAAGGGCT         0     0     0     0     0     0       0       0
D1S_AAACCTGGTAGCGCTC         0     0     0     0     0     0       0       0
                     percent_mt gRNA_umi_count
D1S_AAACCTGAGCGATTCT   3.370787              1
D1S_AAACCTGCAAGCCGTC   5.921053              1
D1S_AAACCTGCAGGGCATA   2.860677              1
D1S_AAACCTGGTAAGGGCT   3.045017              6
D1S_AAACCTGGTAGCGCTC   3.240824              7

TCells_D1N TCells_D1S TCells_D2N TCells_D2S 
      5533       6843       5144       7435 
14278 stimulated cells. 
10677 unstimulated cells.

Run simulated data

dir.create("/project2/xinhe/kevinluo/GSFA/music_analysis/Stimulated_T_Cells/stimulated", recursive = TRUE, showWarnings = FALSE)
res_dir <- "/project2/xinhe/kevinluo/GSFA/music_analysis/Stimulated_T_Cells/stimulated"
# setwd(res_dir)
dir.create(file.path(res_dir,"/music_output"), recursive = TRUE, showWarnings = FALSE)

0. Load input data

feature.names <- data.frame(fread(paste0(data_dir, "GSE119450_RAW/D1N/genes.tsv"),
                                  header = FALSE), stringsAsFactors = FALSE)

expression_profile <- combined_obj@assays$RNA@counts[, stimulated_cells]
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:24]
targets[targets == "NonTarget"] <- "CTRL"
cat("Targets: \n")
print(targets)

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

Dimension of expression profile matrix: 
[1] 33694 14278
Targets: 
 [1] "ARID1A"   "BTLA"     "C10orf54" "CBLB"     "CD3D"     "CD5"     
 [7] "CDKN1B"   "DGKA"     "DGKZ"     "HAVCR2"   "LAG3"     "LCP2"    
[13] "MEF2D"    "CTRL"     "PDCD1"    "RASA2"    "SOCS1"    "STAT6"   
[19] "TCEB2"    "TMEM222"  "TNFRSF9"

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

crop_seq_imputation <- Data_imputation(crop_seq_qc$expression_profile,
                                       crop_seq_qc$perturb_information,
                                       cpu_num = 5)
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 = 5)
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")

crop_seq_vargene <- readRDS("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")

## 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_20_topics.rds")

topic_func <- readRDS("music_output/topic_func_20_topics.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_20_topics.rds")

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

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

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

   perturbation ranking     Score off_target
1          CD3D       1 207.89434       none
2          CBLB       2 181.54109       none
3        HAVCR2       3 164.90364       none
4         MEF2D       4 159.47255       none
5          LAG3       5 146.92054       none
6          LCP2       6 142.39084       none
7       TMEM222       7 131.76594       none
8         RASA2       8 123.32882       none
9         PDCD1       9 117.64996       none
10      TNFRSF9      10 102.44654       none
11          CD5      11 100.44077       none
12        SOCS1      12  97.63884       none
13         BTLA      13  93.85907       none
14         DGKZ      14  84.20973       none
15        TCEB2      15  75.95113       none
16       ARID1A      16  74.05189       none
17       CDKN1B      17  68.22507       none
18        STAT6      18  41.90055       none
19         DGKA      19  39.74996       none
20     C10orf54      20  31.02633       none

calculate the 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_20_topics.rds")

    topic perturbation ranking
1  Topic1         LCP2       1
2  Topic1         CD3D       2
3  Topic1          CD5       3
4  Topic1       ARID1A       4
5  Topic1         DGKA       5
6  Topic1       CDKN1B       6
7  Topic1        STAT6       7
8  Topic1        MEF2D       8
9  Topic1     C10orf54       9
10 Topic1         DGKZ      10

calculate the 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_20_topics.rds")

   Perturbation_1 Perturbation_2 Correlation
2            BTLA         ARID1A  0.81287901
3        C10orf54         ARID1A  0.35807563
23       C10orf54           BTLA  0.16626640
24           CBLB           BTLA -0.53610756
4            CBLB         ARID1A -0.08829456
44           CBLB       C10orf54  0.04880987
65           CD3D           CBLB -0.83212031
25           CD3D           BTLA  0.63583820
5            CD3D         ARID1A  0.18307594
45           CD3D       C10orf54  0.06519028

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

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

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