Master regulator inference

Re-construct TRNs

ARACNE was used to reconstruct transcriptional regulatory netowrks (TRNs) using a public MM gene expression data set GSE6477.

Prepare input files for ARACNE:

library("bcellViper")
data(bcellViper, package="bcellViper")

# load gene expression
load("RData/GSE6477_Symbol.RData")

gene = rownames(ExS)
mT = data.frame(gene, ExS)
# write.table(mT, file = "ARACNE/raw/GSE6477.txt", row.names = FALSE, col.names = TRUE, sep = "\t", quote = F)

# find regulators
Regulator = intersect(names(regulon), gene)
# write.table(cbind(Regulator), file = "ARACNE/raw/Regulator.txt", row.names = FALSE, col.names = FALSE, sep = "\t", quote = F)

Run ARACNE:

Thresholding the mutual information threshold:

java -Xmx8g -jar /sibcb2/bioinformatics/software/ARACNe-AP/Aracne.jar \
    -e ARACNE/raw/GSE6477.txt \
    -o ARACNE/out \
    --tfs ARACNE/raw/Regulator.txt \
    --pvalue 1E-8 \
    --seed 1 \
    --calculateThreshold

Run 100 bootstraps:

for(i in 1:100){

    cat("running bootstrap",i,"\n")
    cmd = paste("java -Xmx8g -jar 
        /sibcb2/bioinformatics/software/ARACNe-AP/Aracne.jar 
        -e ARACNE/raw/GSE6477.txt 
        -o ARACNE/out 
        --tfs ARACNE/raw/Regulator.txt 
        --pvalue 1E-8 
        --seed", i)
    system(cmd)
}

Consolidate all networks:

java -Xmx8g -jar /sibcb2/bioinformatics/software/ARACNe-AP/Aracne.jar \
    -o ARACNE/out/ \
    --consolidate

load ARACNE network

network = read.table(file = "ARACNE/network_consolidated.txt", row.names = NULL, header = TRUE, sep = "\t")
head(network)

  Regulator Target        MI     pvalue
1     FUBP3 TOPORS 0.4333787 0.00000000
2    ZNF593  MRPL4 0.3711633 0.00000000
3     TAF12   MRC1 0.2746556 0.02477282
4      CBX2   PRB3 0.3652976 0.02477282
5      CBX2   PRB4 0.3449795 0.00000000
6      CBX2   PRG3 0.3225342 0.02477282

MM_regul <- aracne2regulon("ARACNE/network_consolidated.txt", ExS)
save(MM_regul, file = "RData/MM_regul.RData")

VIPER

We used VIPER to identify master regulators that potential drive MM progression.

rm(list = ls())
library("ggplot2")
library("viper")
source("/sibcb2/bioinformatics/Script/PathwayEnrichment.R")

# Load Expression data
load("RData/GSE6477_Symbol.RData")
load("RData/MM_regul.RData")
load("RData/rainbowMM_DGE_top.RData")

Enrichment of MM progression signature in patient samples:

# get mean expression
Index = sapply(strsplit(colnames(ExS), "_"), function(x) x[2]) == "NewMM"
Ex    = ExS[, Index]
refMean    <- apply(Ex, 1, mean)
ExCentered <- Ex - refMean

# Load expression signature
UP <- sigList$MM1S$UP
DN <- sigList$MM1S$DN

# get Z-score for the enrichment of signatures
M  <- ncol(ExCentered)
zPOS <- rep(0, M)
zNEG <- rep(0, M)

for(i in 1:M){
  
  fcVector <- ExCentered[,i]
  gsList   <- list(UP = UP, DN = DN)
  
  tempZ    <- SimpleRankTest(fcVector, gsList)
  zPOS[i]  <- tempZ[1]
  zNEG[i]  <- tempZ[2]
}
zDIF        <- zPOS - zNEG
names(zDIF) <- colnames(Ex)
zDIF        <- sort(zDIF)

Comparing samples that are more agressive and those that are less agressive:

# get samples with high or low z-scores
cutoff = median(zDIF)
hName  = names(zDIF)[zDIF > cutoff]
lName  = names(zDIF)[zDIF < cutoff]
hM     = ExS[, hName]
lM     = ExS[, lName]
group  = c(rep("high", length(hName)), rep("low", length(lName)))
signature = rowTtest(hM, lM)
signature <- (qnorm(signature$p.value/2, lower.tail = FALSE) *sign(signature$statistic))[, 1]
nullmodel = ttestNull(hM, lM, per = 1000, repos = TRUE, seed = 1, verbose = FALSE)
mrs <- msviper(signature, MM_regul, nullmodel, verbose = FALSE)

Plot top master regulators:

Tx = summary(mrs, 100)
# write.table(Tx, file = "ARACNE/Enriched_regulators.txt", row.names = FALSE, col.names = TRUE, sep = "\t", quote = FALSE)
options(repr.plot.width=6, repr.plot.height=10, repr.plot.res = 200)
plot(mrs, 30, cex = 1)