seeds-reproducibility RWR

examples/DTSEA_comparison/DTSEA_test.R

@@ -14,7 +14,7 @@ if (!require(pacman, quietly = TRUE)){
 
 
 library(pacman)
-p_load(magrittr, dplyr, BiocManager, devtools, here, ggplot2)
+p_load(magrittr, dplyr, BiocManager, devtools, here, ggplot2, tibble)
 
 
 # Install the required packages for DTSEA and DisGeNet
@@ -51,29 +51,31 @@ colnames(drug_targets)<- gsub("drugbank_id", "drug_id", colnames(drug_targets))
 ##############
 # Perform a simple DTSEA analysis using default optional parameters then sort
 # the result dataframe by normalized enrichment scores (NES)
+set.seed(234)
 result_FA <- DTSEA(network = graph,
                 disease = FA_signature$Gene,
                 drugs = drug_targets, verbose = FALSE) %>% arrange(desc(NES))
 
 relevantFA_results <- select(result_FA, -leadingEdge) %>% arrange(desc(NES)) %>%  filter(NES > 0 & padj < .01)
-NES_FAdrug_targets <- relevantFA_results[,c(1,2,6)] %>%  add_column(target = drug_targets$gene_target[match(relevantFA_results$drug_id, drug_targets$drug_id)]) %>%  
+NES_FAdrug_targets <- relevantFA_results[,c(1,3,6)] %>%  add_column(target = drug_targets$gene_target[match(relevantFA_results$drug_id, drug_targets$drug_id)]) %>%  
   add_column(drug = drug_targets$drug_name[match(relevantFA_results$drug_id, drug_targets$drug_id)]) 
-length(unique(NES_FAdrug_targets$target)) ## 54 unique targets
+length(unique(NES_FAdrug_targets$target)) ## 52 unique targets
 
 
 ############
 ## DTSEA: FM use-case
 ##############
 # Perform a simple DTSEA analysis using default optional parameters then sort
 # the result dataframe by normalized enrichment scores (NES)
+set.seed(234)
 result_FM <- DTSEA(network = graph,
                    disease = FM_signature$symbol,
                    drugs = drug_targets, verbose = FALSE) %>% arrange(desc(NES))
 
 relevantFM_results <- select(result_FM, -leadingEdge) %>% arrange(desc(NES)) %>%  filter(NES > 0 & padj < .01)
-NES_FMdrug_targets <- relevantFM_results[,c(1,2,6)] %>%  add_column(target = drug_targets$gene_target[match(relevantFM_results$drug_id, drug_targets$drug_id)]) %>%  
+NES_FMdrug_targets <- relevantFM_results[,c(1,3,6)] %>%  add_column(target = drug_targets$gene_target[match(relevantFM_results$drug_id, drug_targets$drug_id)]) %>%  
   add_column(drug = drug_targets$drug_name[match(relevantFM_results$drug_id, drug_targets$drug_id)])
-length(unique(NES_FMdrug_targets$target)) ## 85 unique targets
+length(unique(NES_FMdrug_targets$target)) ## 91 unique targets
 
 
 ###############################
@@ -87,12 +89,12 @@ length(unique(NES_FMdrug_targets$target)) ## 85 unique targets
 ### 1. Load relevant filtered drexml results 
 FA_drexml<- read.delim(gzfile(here("examples", "fanconi_anemia" , "results","shap_filtered_stability_symbol.tsv.gz"))) %>% column_to_rownames("circuit_name")
 
-overlapFA <- colnames(FA_drexml)[colnames(FA_drexml) %in% NES_FAdrug_targets$target] ## 7 overlap
+overlapFA <- colnames(FA_drexml)[colnames(FA_drexml) %in% NES_FAdrug_targets$target] ## 6 overlap
 
 ### 2. Rank targets in NES_FAdrug_targets based on p-value
 nes_fa_ranks <- NES_FAdrug_targets  %>%
   mutate(rank_nes_fa = row_number()) %>%
-  select(drug_id, drug ,target, pval, NES, rank_nes_fa) %>%
+  select(drug_id, drug ,target, padj, NES, rank_nes_fa) %>%
   filter(target %in% overlapFA)
 
 ### 3. Calculate mean absolute SHAP values and count relevant circuits for all targets
@@ -109,7 +111,7 @@ FA_stats_df <- as.data.frame(t(FA_stats), stringsAsFactors = FALSE) %>%
 
 ### 4. Create summarisation table
 
-summary_DTSEAcomparisonFA <- merge(nes_fa_ranks, FA_stats_df, by = "target", all.x = TRUE)  %>% .[order(.$pval,decreasing = F),]
+summary_DTSEAcomparisonFA <- merge(nes_fa_ranks, FA_stats_df, by = "target", all.x = TRUE)  %>% .[order(.$padj,decreasing = F),]
 write.table( summary_DTSEAcomparisonFA, here("examples", "DTSEA_comparison" ,"summary_DTSEAcomparisonFA.tsv"), sep="\t", quote = F, row.names = F, col.names = T)
 
 ### 5. Compare with cmap results
@@ -130,7 +132,7 @@ overlapFM <- colnames(FM_drexml)[colnames(FM_drexml) %in% NES_FMdrug_targets$tar
 ### 2. Rank targets in NES_FAdrug_targets based on p-value
 nes_fm_ranks <- NES_FMdrug_targets  %>%
   mutate(rank_nes_fm = row_number()) %>%
-  select(drug_id, drug ,target, pval, NES, rank_nes_fm) %>%
+  select(drug_id, drug ,target, padj, NES, rank_nes_fm) %>%
   filter(target %in% overlapFM)
 
 ### 3. Calculate mean absolute SHAP values and count relevant circuits for all targets
@@ -147,7 +149,7 @@ FM_stats_df <- as.data.frame(t(FM_stats), stringsAsFactors = FALSE) %>%
 
 ### 4. Create summarisation table
 
-summary_DTSEAcomparisonFM <- merge(nes_fm_ranks, FM_stats_df, by = "target", all.x = TRUE) %>% .[order(.$pval,decreasing = F),]
+summary_DTSEAcomparisonFM <- merge(nes_fm_ranks, FM_stats_df, by = "target", all.x = TRUE) %>% .[order(.$padj,decreasing = F),]
 write.table( summary_DTSEAcomparisonFM, here("examples", "DTSEA_comparison" ,"summary_DTSEAcomparisonFM.tsv"), sep="\t" ,quote = F, row.names = F, col.names = T)
 
 ### 5. Compare with cmap results

examples/DTSEA_comparison/summary_DTSEAcomparisonFA_curated.csv

@@ -0,0 +1,9 @@
+target	drug_id	drug	padj	NES	rank_nes_overall	mean_abs_shap	circuits_above_zero
+PIK3CD	DB11891	Fimepinostat	5.86409023807469E-05	2.35164642871307	17	0.256086674697682	88
+EGFR	DB12174	CUDC-101	0.000107047195274	2.30956196828472	18	0.00192081948265	4
+VEGFA	DB10772	Foreskin keratinocyte (neonatal)	0.000605901166843	2.25022081918057	22	0.00120908257817	1
+TNF	DB05992	Plinabulin	0.001806398537351	2.17902906869224	34	0.92890530672613	117
+EGFR	DB03496	Alvocidib	0.001922283066699	2.1385974800089	43	0.00192081948265	4
+LPL	DB13751	Glycyrrhizic acid	0.002994358115833	1.98861632676143	80	0.119130985152545	41
+TUBB3	DB04845	Ixabepilone	0.003973089808473	2.12724378082689	46	0.059374744208029	19
+VEGFA	DB05294	Vandetanib	0.006768006463629	2.06043311812418	65	0.00120908257817	1

examples/DTSEA_comparison/summary_DTSEAcomparisonFM_curated.csv

@@ -0,0 +1,17 @@
+target	drug_id	drug	padj	NES	rank_nes_overall	mean_abs_shap	circuits_above_zero
+LDHB	DB11638	Artenimol	1.98408865659993E-07	2.40200990483904	28	0.053377641079268	5
+RRM1	DB00631	Clofarabine	2.69329582279362E-07	2.49684342163065	20	0.006174314357354	1
+EGFR	DB03496	Alvocidib	7.07938283343453E-07	2.41138188947284	26	0.169799002032272	14
+PSMB2	DB11762	Marizomib	1.60770306008383E-06	2.49958816383425	19	0.005816669406011	1
+PDGFRB	DB09283	Trapidil	3.65221302466481E-06	2.41975957955149	24	0.023702379013807	3
+AR	DB00421	Spironolactone	0.000103160770396	2.25020912161814	66	0.005380756008407	1
+EGFR	DB12174	CUDC-101	0.000123161358813	2.20658887402019	72	0.169799002032272	14
+KDR	DB05608	MKC-1	0.00102343053143	2.1573956192271	80	0.00652201980873	1
+PDGFRB	DB12147	Erdafitinib	0.001183330282605	2.09775818402175	94	0.023702379013807	3
+CACNA2D1	DB00230	Pregabalin	0.001229636512096	2.09132430764859	97	0.100578160553848	9
+CACNA2D1	DB08872	Gabapentin enacarbil	0.001229636512096	2.09132430764859	98	0.100578160553848	9
+AR	DB08804	Nandrolone decanoate	0.00195639906786	2.08835351757494	104	0.005380756008407	1
+PDGFRB	DB11694	Ilorasertib	0.003464725355424	2.00282691528826	135	0.023702379013807	3
+AR	DB00717	Norethisterone	0.007142878402106	1.66398775624387	238	0.005380756008407	1
+AR	DB08867	Ulipristal	0.007142878402106	1.66398775624387	239	0.005380756008407	1
+AR	DB14583	Segesterone acetate	0.007142878402106	1.66398775624387	240	0.005380756008407	1