update hardcoded 135

sceptre_paper/analysis_drivers/analysis_drivers_xie/cis_enrichment_analysis_9.R

@@ -45,6 +45,7 @@ resampling_results = resampling_results_xie_cis
 ss_xie_cis = readRDS(file = paste0(processed_dir, '/ss_xie_cis.rds'))
 original_results <- ss_xie_cis %>% select('gene_id', 'gRNA_id', 'ss.down', 'reject.down') %>% dplyr::rename(rejected = reject.down)
 
+# BUG! ROWS NOT ORDERED IDENTICALLY.
 original_results = cbind(original_results, resampling_results[, c('gene_short_name', 'chr', 'strand', 'target_gene.start',
                                                                   'target_gene.stop', 'TSS', 'target_site.start', 'target_site.stop',
                                                                   'target_site.mid')])

sceptre_paper/analysis_drivers/analysis_drivers_xie/pre_process_data_3.R

@@ -28,6 +28,8 @@ h5_odm <- create_ondisc_matrix_from_h5_list(h5_list = h5_files,
                                             metadata_fp = paste0(processed_dir, "/odm/gene_expression_metadata.rds"),
                                             barcode_suffixes = batch_lane,
                                             progress = TRUE)
+h5_odm <- read_odm(odm_fp = paste0(processed_dir, "/odm/gene_expression_matrix.odm"),
+                   metadata_fp = paste0(processed_dir, "/odm/gene_expression_metadata.rds"))
 
 # modify the ondisc matrix by identifying the protein-coding genes
 gene_df <- read.xlsx(xlsxFile = paste0(raw_data_dir, "/Genes.xlsx"), sheet = 1)

sceptre_paper/analysis_drivers/analysis_drivers_xie/significance_score_8.R

@@ -83,6 +83,10 @@ ss_xie_cis = data.frame(gene_id = gRNA.gene.pair$gene_id[gRNA.gene.pair$type ==
                         gRNA_id = gRNA.gene.pair$gRNA_id[gRNA.gene.pair$type == 'cis'],
                         ss.down = SS.down.cis[cbind(m.gene.id, m.gRNA.id)])
 
-ss_thres = sort(ss_xie_cis$ss.down, decreasing = T)[135]
+# load the number of cis SCEPTRE discoveries
+sceptre_res <- fst::read_fst(paste0(results_dir, "/all_results_annotated.fst"))
+n_sceptre_discoveries <- sceptre_res %>% dplyr::filter(type == "cis") %>% dplyr::pull(rejected) %>% sum()
+
+ss_thres = sort(ss_xie_cis$ss.down, decreasing = T)[n_sceptre_discoveries]
 ss_xie_cis$reject.down = ss_xie_cis$ss.down >= ss_thres
 saveRDS(ss_xie_cis, file = paste0(processed_dir, '/ss_xie_cis.rds'))

sceptre_paper/plotting/Figure5.R

@@ -1,5 +1,6 @@
 args <- commandArgs(trailingOnly = TRUE)
-code_dir <- if (is.na(args[1])) "/Users/timbarry/Box/SCEPTRE/SCEPTRE/" else args[1]
+code_dir <- paste0(.get_config_path("LOCAL_CODE_DIR"), "sceptre-manuscript")
+
 library(mgcv)
 library(scales)
 require(katsevich2020)
@@ -28,7 +29,7 @@ table(rejected_by_annotated)
 #   Neither method      Both methods      SCEPTRE only Virtual FACS only 
 #             3367                84                51                51 
 
-ss_thres =sort(ss_xie_cis$ss.down, decreasing = T)[135]
+ss_thres = sort(ss_xie_cis$ss.down, decreasing = T)[135]
 df_s = data.frame(gene_id = resampling_results$gene_id, gRNA_id = resampling_results$gRNA_id, 
                   p_value = resampling_results$p_value, signif.score = ss_xie_cis$ss.down, 
                   rejected_by_annotated = rejected_by_annotated)

sceptre_paper/plotting/load_data_for_plotting.R

@@ -1,4 +1,5 @@
-offsite_dir <- if (is.na(args[2])) "/Volumes/tims_new_drive/research/sceptre_files/" else args[2]
+# offsite_dir <- if (is.na(args[2])) "/Volumes/tims_new_drive/research/sceptre_files/" else args[2]
+offsite_dir <- .get_config_path("LOCAL_SCEPTRE_DATA_DIR")
 manuscript_figure_dir <- paste0(code_dir, "/sceptre_paper/manuscript/figures")
 
 # Gasperini results
@@ -10,14 +11,17 @@ covariates_gasp = read.fst(sprintf("%s/data/gasperini/processed/cell_covariate_m
 # Xie results: 4 methods: original, resampling_results (sceptre), likelihood results (iNB), and monocle_nb
 gRNA.gene.pair = read.fst(sprintf("%s/data/xie/processed/gRNA_gene_pairs.fst", offsite_dir)) %>% as_tibble()
 original_results_xie = readRDS(sprintf("%s/data/xie/processed/raw_pval_xie.rds", offsite_dir)) %>% as_tibble()
+# original_results_xie = readRDS("~/Desktop/raw_pval_xie.rds") %>% as_tibble()
 resampling_results_xie_with_names = read.fst(sprintf("%s/results/xie/sceptre/all_results_annotated.fst", offsite_dir)) %>% as_tibble()
+# resampling_results_xie_with_names = read.fst("~/Desktop/all_results_annotated.fst") %>% as_tibble()
 likelihood_results_xie = read.fst(sprintf("%s/results/xie/negative_binomial/all_results.fst", offsite_dir)) %>% as_tibble()
 
 p_vals_bulk <- paste0(offsite_dir, "/results/xie/bulk_rna_seq/pvals_arl15_enh.rds") %>% readRDS() %>% as_tibble() %>% rename(gene_names = gene_id)
 p_vals_bulk_myb3_enh3 <- paste0(offsite_dir, "/results/xie/bulk_rna_seq/pvals_myb_enh3.rds") %>% readRDS() %>% as_tibble() %>% rename(gene_names = gene_id)
 covariates_xie = read.fst(sprintf("%s/data/xie/processed/covariate_model_matrix.fst", offsite_dir)) %>% as_tibble()
 grna_indicator_matrix_xie = read.fst(sprintf("%s/data/xie/processed/gRNA_indicator_matrix.fst", offsite_dir)) %>% as_tibble()
 ss_xie_cis = readRDS(sprintf("%s/data/xie/processed/ss_xie_cis.rds", offsite_dir)) %>% as_tibble()
+# ss_xie_cis = readRDS("~/Desktop/ss_xie_cis.rds") %>% as_tibble()
 
 resampling_results_xie_cis = resampling_results_xie_with_names %>% filter(type == 'cis')
 gene.mart = readRDS(sprintf("%s/data/xie/processed/gene_mart.rds", offsite_dir))