Loading MRs

Master regulators were predicted from multiple AML data sets, as it was shown for TCGA data.

library("dplyr")
library("ggplot2")
library("pheatmap")

mergedT = read.table(file = "/sibcb2/bioinformatics/LSC/out/AML_MR_17gene_summary.txt", 
  row.names = 1, header = TRUE, sep = "\t")

MR = as.character(read.table(file = "/sibcb2/bioinformatics/LSC/out/AML_MR_17gene_PG.txt", 
  row.names = NULL, header = FALSE, sep = "\t")[,1])

cMatrix = matrix(0, length(MR), 2)
dimnames(cMatrix) = list(MR, c("Activated", "Repressed"))
cMatrix[rownames(mergedT), ] = as.matrix(mergedT)

MR_UP  <- rownames(filter(mergedT, Activated > 0))
NumUP  <- cMatrix[MR, "Activated"]

MR_DN  <- rownames(filter(mergedT, Repressed > 0))
NumDN  <- cMatrix[MR, "Repressed"]

CRISPR hits

load("DataSummary/Focused_Depleted.RData")

MA_vitro  <- maList_HS$vitro_vs_plasmid
MA_vivo   <- maList_HS$Post_vs_Pre
HM_vitro  <- hmList_HS$vitro_vs_plasmid
HM_vivo   <- hmList_HS$Post_vs_Pre
CRISPR    <- list(MA_vitro = MA_vitro, MA_vivo = MA_vivo, HM_vitro = HM_vitro, HM_vivo = HM_vivo)

Index_maVitro = (MR %in% MA_vitro) + 0
Index_maVivo  = (MR %in% MA_vivo)  + 0
Index_hmVitro = (MR %in% HM_vitro) + 0
Index_hmVivo  = (MR %in% HM_vivo)  + 0
Index_Focused = (MR %in% PG_HS) + 0
mergedTable = data.frame(MR, NumUP, NumDN, Index_maVitro, Index_maVivo, Index_hmVitro, Index_hmVivo, Index_Focused)
#write.table(mergedTable, file = "out/VIPER_summary_Focused.txt", row.names = FALSE, col.names = TRUE, sep = "\t", quote = FALSE)

Enrichment

res = matrix(NA, 4, 6)
PG      = intersect(MR, PG_HS)
for(i in 1:4){

  temp = intersect(CRISPR[[i]], PG)
  NC   = setdiff(PG, temp)

  a = length(intersect(MR_UP, temp))
  b = length(intersect(MR_UP, NC))
  c = length(intersect(MR_DN, temp))
  d = length(intersect(MR_DN, NC))
  abcd = matrix(c(a,b,c,d), 2,2, byrow = TRUE)
  dimnames(abcd) = list(c("MR_UP", "MR_DN"), c("Depleted", "NC"))

  res[i, 1:4] = c(a, b, c, d)
  res[i, 5  ] = fisher.test(abcd)$p.value
  res[i, 6  ] = fisher.test(abcd)$estimate
}
dimnames(res) = list(names(CRISPR), c("UP_Depleted", "UP_NotDepleted", "DN_Depleted", "DN_NotDepleted", "Pvalue", "OR"))
knitr::kable(data.frame(res))
UP_Depleted UP_NotDepleted DN_Depleted DN_NotDepleted Pvalue OR
MA_vitro 34 18 12 18 0.0373880 2.795665
MA_vivo 31 21 16 14 0.6462778 1.287595
HM_vitro 35 17 12 18 0.0210834 3.043321
HM_vivo 36 16 17 13 0.3379715 1.708946 
res = dplyr::mutate(data.frame(res), 
  Fra_UP = UP_Depleted/(UP_Depleted + UP_NotDepleted),
  Fra_DN = DN_Depleted/(DN_Depleted + DN_NotDepleted))

out = data.frame(name= rep(rownames(res), 2), 
  Fra = c(res$Fra_UP, res$Fra_DN), 
  group = rep(c("UP", "DN"), each = 4))

  p <- ggplot(out, aes(x = name, y = Fra, fill = group)) + ylim(0, 1)
  p <- p + theme_bw() + labs(x = "CRISPR conditions", y = "Depleted fraction", title = "") 
  p <- p + geom_bar(stat="identity", position="dodge")
  p <- p + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
  p <- p + scale_fill_manual(values = c("UP" = "#dd1c77", "DN" = "grey"))
  print(p)

Activated hits

out = filter(mergedTable, MR %in% PG_HS, NumUP >= 1)
s1 = out[, "Index_maVitro"] + out[, "Index_hmVitro"]
s2 = out[, "Index_maVivo"]  + out[, "Index_hmVivo"]
group = rep("NC", nrow(out))
group[s1 >  0 & s2 >  0] = "shared"
group[s1 >  0 & s2 == 0] = "vitro"
group[s1 == 0 & s2 >  0] = "vivo"
out = data.frame(out, group)
t1 = filter(out, group == "shared") %>% arrange(desc(Index_maVitro + Index_hmVitro + Index_maVivo + Index_hmVivo))
t2 = filter(out, group == "vitro") %>% arrange(desc(Index_maVitro + Index_hmVitro + Index_maVivo + Index_hmVivo))
t3 = filter(out, group == "vivo") %>% arrange(desc(Index_maVitro + Index_hmVitro + Index_maVivo + Index_hmVivo))
t4 = filter(out, group == "NC") %>% arrange(desc(Index_maVitro + Index_hmVitro + Index_maVivo + Index_hmVivo))
oTable = rbind(t1, t2, t3, t4)
reM = as.matrix(oTable[,c("Index_maVitro", "Index_maVivo", "Index_hmVitro", "Index_hmVivo")])
colnames(reM) = gsub("V", "_V", sapply(strsplit(colnames(reM), "_"), function(x) x[2]))

annotation_row = oTable[, "group", drop = FALSE]
ann_colors = list(reGroup = c(shared = "#7b3294", vitro = "#c2a5cf", vivo = "#a6dba0", NC = "#008837"))

pheatmap(reM,
         main = "",
         annotation_row = annotation_row,
         annotation_colors = ann_colors,         
         border_color = "grey60", color = c("#f5f5f5","#f5f5f5","#dd1c77"),
         legend_breaks = c(-1,0,1),
         legend_labels = c("NC", "NC", "Depleted"),
         fontsize_row = 6,
         show_rownames = TRUE, show_colnames = TRUE,
         cluster_rows  = FALSE, cluster_cols  = FALSE)