In this document, we will show how to predict master regulators that potential drive AML fitness using TCGA AML gene expression data.

Patient stratification

library("ggplot2")
library("viper")
source("/sibcb2/bioinformatics/Script/PathwayEnrichment.R")

load("/sibcb2/bioinformatics/KnowledgeBase/Firehose_Expression/Firehose_Expression_Tumor.RData")
load("/sibcb2/bioinformatics/Placenta/RData/TCGA_regulonDB.RData")
load("/sibcb2/bioinformatics/LSC/RData/Signature_Top500.RData")

pw      = list(up = Signature_Top500$sigUP_Hs, dn = Signature_Top500$sigDN_Hs)
regulon = TCGA_regulonDB[["regulonlaml"]]
Ex      = Firehose_Expression_Tumor[["LAML"]]
ref     = apply(Ex, 1, mean)
dM      = Ex - ref

zscore  = rep(NA, ncol(dM))
for(j in 1:ncol(dM)){

  temp = SimpleRankTest(dM[,j], pw)
  zscore[j] = temp[1] - temp[2]
}
names(zscore) = colnames(dM)

Enrichment = rep("others", length(zscore))
Enrichment[zscore >= quantile(zscore, 1 - 50/length(zscore))] = "high"
Enrichment[zscore <= quantile(zscore,     50/length(zscore))] = "low"
temp = data.frame(rnk = rank(zscore), zscore, Enrichment)

p <- ggplot(temp, aes(x = rnk, y = zscore, color = Enrichment, fill = Enrichment)) + geom_bar(stat="identity") + theme_bw()
p <- p + theme(plot.title = element_text(hjust = 0.5))
p <- p + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + ylab("Enrichment z-score") + xlab("") + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())
p <- p + scale_colour_manual(values = c("red", "blue", "grey"))
print(p)

MR prediction

ExH = Ex[, names(sort(zscore, decreasing = TRUE)[1:50])]
ExL = Ex[, names(sort(zscore, decreasing = FALSE )[1:50])]

signature = rowTtest(ExH, ExL)
signature = (qnorm(signature$p.value / 2, lower.tail = FALSE) * sign(signature$statistic))[, 1]
nullmodel = ttestNull(ExH, ExL, per = 500, repos = TRUE, verbose = FALSE)
mrs     = msviper(signature, regulon, nullmodel, verbose = FALSE)
knitr::kable(data.frame(summary(mrs, 10)))
Regulon Size NES p.value FDR
YWHAE YWHAE 95 4.91 9.00e-07 0.00273
SDHD SDHD 89 4.65 3.30e-06 0.00411
ZBTB8OS ZBTB8OS 94 4.60 4.20e-06 0.00411
PRDM15 PRDM15 72 -4.29 1.76e-05 0.01040
LCK LCK 33 -4.45 8.80e-06 0.00574
SEPT5 SEPT5 108 -4.46 8.10e-06 0.00574
EPHA4 EPHA4 45 -4.48 7.40e-06 0.00574
KLF8 KLF8 137 -4.60 4.10e-06 0.00411
RUNDC3A RUNDC3A 70 -4.79 1.60e-06 0.00325
TREML1 TREML1 67 -4.92 9.00e-07 0.00273

sessionInfo

sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-conda_cos6-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server 7.7 (Maipo)

Matrix products: default
BLAS: /sibcb2/bioinformatics/software/BcbioNG/anaconda/lib/R/lib/libRblas.so
LAPACK: /sibcb2/bioinformatics/software/BcbioNG/anaconda/lib/R/lib/libRlapack.so

locale:
[1] C

attached base packages:
[1] parallel  grid      stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] viper_1.16.0        scales_1.1.0        pheatmap_1.0.12    
 [4] dplyr_1.0.5         survcomp_1.32.0     prodlim_2019.11.13 
 [7] RColorBrewer_1.1-2  gridExtra_2.3       survHD_0.99.1      
[10] survC1_1.0-2        Hmisc_4.3-0         Formula_1.2-3      
[13] lattice_0.20-38     penalized_0.9-51    survival_3.1-8     
[16] Biobase_2.42.0      BiocGenerics_0.28.0 VennDiagram_1.6.20 
[19] futile.logger_1.4.3 ggrepel_0.8.1       ggplot2_3.2.1      
[22] rmarkdown_2.0      

loaded via a namespace (and not attached):
 [1] mixtools_1.1.0       splines_3.5.1        rmeta_3.0           
 [4] highr_0.8            latticeExtra_0.6-28  bootstrap_2019.6    
 [7] yaml_2.2.0           survivalROC_1.0.3    pillar_1.5.1        
[10] backports_1.1.5      glue_1.4.2           digest_0.6.23       
[13] checkmate_1.9.4      colorspace_1.4-1     htmltools_0.5.1.1   
[16] Matrix_1.2-18        pkgconfig_2.0.3      purrr_0.3.3         
[19] lava_1.6.7           htmlTable_1.13.3     tibble_3.1.0        
[22] generics_0.0.2       farver_2.0.1         ellipsis_0.3.0      
[25] withr_2.1.2          nnet_7.3-12          lazyeval_0.2.2      
[28] magrittr_1.5         crayon_1.3.4         evaluate_0.14       
[31] fansi_0.4.0          segmented_1.1-0      MASS_7.3-51.5       
[34] class_7.3-15         foreign_0.8-74       SuppDists_1.1-9.5   
[37] tools_3.5.1          data.table_1.12.8    formatR_1.7         
[40] lifecycle_1.0.0      stringr_1.4.0        munsell_0.5.0       
[43] cluster_2.1.0        lambda.r_1.2.4       e1071_1.7-3         
[46] compiler_3.5.1       rlang_0.4.10         rstudioapi_0.10     
[49] htmlwidgets_1.5.1    base64enc_0.1-3      labeling_0.3        
[52] gtable_0.3.0         DBI_1.1.0            R6_2.4.1            
[55] knitr_1.26           utf8_1.1.4           futile.options_1.0.1
[58] KernSmooth_2.23-16   stringi_1.4.3        Rcpp_1.0.3          
[61] vctrs_0.3.7          rpart_4.1-15         acepack_1.4.1       
[64] tidyselect_1.1.0     xfun_0.18