Code to reproduce all figures

library(iTOP)
library(pcalg)
library(readr)
library(bcv)
library(NMF)
library(igraph)
library(Rgraphviz)
library(foreach)
library(doMC)
library(glmnet)
library(TANDEM)
library(latex2exp)
library(beeswarm)

If you have trouble installing pcalg, please use the following code (as sometimes installing dependencies from Bioconductor seems to fail for this package):

dependencies = c("Rgraphviz", "graph", "RBGL")
source("https://bioconductor.org/biocLite.R")
for(dependency in dependencies) {
  if(!require(dependency, character.only=T)) {
    biocLite(dependency)
  }
}
install.packages("pcalg")

In addition, on some Linux systems, we also needed to install some additional libraries:

• sudo apt-get install libv8-3.14-dev
• sudo apt-get install libgmp3-dev

Figure 1

Is actually a cartoon that is not based on actual data, but rather meant as a high-level description of the manuscript. Hence, there is no code for figure 1.

Figure 2

A very simple example of using the RV coefficient to compare datasets.

col = rep(c("#ca5670","#72a555"), each=50)

set.seed(1)
n = 100
p = 2
x = matrix(rnorm(n*p),n,p)
x[1:50,] = x[1:50,] - 10
y = matrix(rnorm(n*p),n,p)
y[1:50,1] = y[1:50,1] + 10
y[1:50,2] = y[1:50,2] - 10
z = matrix(rnorm(n*p),n,p)

x = scale(x, scale=F)
y = scale(y, scale=F)
z = scale(z, scale=F)

data = list(x, y, z)
xlim = c(min(x[,1], y[,1], z[,1]), max(x[,1], y[,1], z[,1]))
ylim = c(min(x[,2], y[,2], z[,2]), max(x[,2], y[,2], z[,2]))
q = length(data)
for(i in 1:q) {
  plot(data[[i]], xlab="", ylab="", ylim=ylim, xlim=xlim, col=col, pch=16, xaxt="n", yaxt="n")
}

config_matrices = compute.config.matrices(data, mod.rv=F)

ann = data.frame(Var1=factor(col, levels=unique(col)))
annColors = list(Var1=unique(col))
for(i in 1:q) {
  aheatmap(config_matrices[[i]], Rowv=NA, Colv=NA, cexRow=0, cexCol=0, breaks=0, annCol=ann, annRow=ann, annLegend=F, annColors=annColors)
}

Figure 3

A simple example of a partial matrix correlation. The last figure shows that the inferred structure (using the PC algorithm) indeed recovers the structure underlying the data (x -> y -> z).

col = rep(c("#ca5670","#72a555","#ab62c0","#c57c3c","#638ccc"), 10)

set.seed(1)
create.data = function(x) {
  n = nrow(x)
  p = ncol(x)
  x2 = NULL
  for(i in 1:10) {
    x2 = rbind(x2, x+0.1*matrix(rnorm(n*p),n,p))
  }
  return(x2)
}

n = 5
p = 2
x = create.data(t(matrix(c(-1, 1, 1, 1, -1, -1, 0, 2, 1, -1), p, n)))
y = create.data(t(matrix(c(-1, 1, 1, 1, 0, -1, 0, 2, 1, -1), p, n)))
z = create.data(t(matrix(c(-1, 1, 1, 1, 0, -1, 0, 2, -1, 0), p, n)))

data = list(x, y, z)
xlim = c(min(x[,1], y[,1], z[,1]), max(x[,1], y[,1], z[,1]))
ylim = c(min(x[,2], y[,2], z[,2]), max(x[,2], y[,2], z[,2]))
q = length(data)
for(i in 1:q) {
  plot(data[[i]], xlab="", ylab="", ylim=ylim, xlim=xlim, col=col, pch=16)
}

config_matrices = compute.config.matrices(data)
cors = rv.cor.matrix(config_matrices)
nperm = nboots = 1000
cors_boot = run.bootstraps(config_matrices, nboots)
cors_perm = run.permutations(config_matrices, nperm)
suffStat = list(cors=cors, cors_perm=cors_perm)

fit = pc(suffStat=suffStat, indepTest=rv.link.significance, labels=c("x", "y", "z"), alpha=0.01, verbose=F, solve.confl=T, conservative=T)
plot(fit, main="")

rv.pcor(cors, 1, 3, 2)

## [1] 0.005361244

rv.conf.interval(cors_boot, 1, 3, 2, .99)

##       0.5%      99.5% 
## -0.2660228  0.2846492

rv.pval(cors, cors_perm, 1, 3, 2)

## [1] 0.354

rv.pcor(cors, 1, 2, 3)

## [1] 0.8777891

rv.conf.interval(cors_boot, 1, 2, 3, .99)

##      0.5%     99.5% 
## 0.7705777 0.9422735

rv.pval(cors, cors_perm, 1, 2, 3)

## [1] 0

rv.pcor(cors, 2, 3, 1)

## [1] 0.3695893

rv.conf.interval(cors_boot, 2, 3, 1, .99)

##      0.5%     99.5% 
## 0.2054663 0.5217131

rv.pval(cors, cors_perm, 2, 3, 1)

## [1] 0

Figure 3

Artificial data experiments that demonstrate the behaviour of the RV coefficient when comparing:

• Two continuous data blocks
• Two binary data blocks

Both with and without centering. Of course, for binary data blocks, Jaccard similarity was used to create the configuration matrices.

rv.coef.wrapper = function(x1, x2, similarity_fun_x1, similarity_fun_x2, center) {
  S1 = compute.config.matrix(x1, similarity_fun_x1, center)
  S2 = compute.config.matrix(x2, similarity_fun_x2, center)
  rv = rv.coef(S1, S2)
  return(rv)
}

cont.cont = function() {
  n = 100
  p = 100
  set.seed(1)
  x = matrix(rnorm(n*p), n, p) + 10
  y = matrix(rnorm(n*p), n, p)
  
  alphas = seq(0, 1, length.out=11)
  res = matrix(NA, length(alphas), 2)
  for(i in 1:length(alphas)) {
    alpha = alphas[i]
    z = (1-alpha)*x + alpha*y
    res[i,1] = rv.coef.wrapper(x, z, inner.product, inner.product, center=T)
    res[i,2] = rv.coef.wrapper(x, z, inner.product, inner.product, center=F)
  }
  plot(alphas, res[,1], ylim=c(min(res),max(res)), ylab="", xlab="", pch=16)
  points(alphas, res[,2], col=2, pch=16)
}

binary.binary = function() {
  n = 100
  p = 100
  set.seed(1)
  x = matrix(rbinom(n*p, 1, 0.5), n, p)
  
  alphas = seq(0, 0.5, length.out=11)
  res = matrix(NA, length(alphas), 2)
  for(i in 1:length(alphas)) {
    alpha = alphas[i]
    y = x
    ind = sample(1:length(y), alpha*length(y), replace=F)
    y[ind] = -y[ind]+1
    res[i,1] = rv.coef.wrapper(x, y, jaccard, jaccard, center=T)
    res[i,2] = rv.coef.wrapper(x, y, jaccard, jaccard, center=F)
  }
  ylim = c(min(res),max(res))
  plot(alphas, res[,1], ylim=ylim, ylab="", xlab="", pch=16)
  points(alphas, res[,2], col=2, pch=16)
}

cont.cont()

binary.binary()

Figure 4

Parsing

The pharmacogenomics dataset analysis. First, let’s parse the data. The required files can be downloaded from:

• Mutations: http://www.cancerrxgene.org/gdsc1000/GDSC1000_WebResources//Data/BEMs/CellLines/CellLines_CG_BEMs.zip (unzip and copy PANCAN_SEQ_BEM.txt)
• CNAs: http://www.cancerrxgene.org/gdsc1000/GDSC1000_WebResources//Data/BEMs/CellLines/CellLines_CNV_BEMs.zip (unzip and copy PANCAN_CNA_BEM.rdata.txt)
• Methylation: http://www.cancerrxgene.org/gdsc1000/GDSC1000_WebResources//Data/preprocessed/methylation/METH_CELL_DATA.txt.zip (unzip and copy F2_METH_CELL_DATA.txt)
• Methylation annotation: http://www.cancerrxgene.org/gdsc1000/GDSC1000_WebResources//Data/otherAnnotations/methSampleId_2_cosmicIds.xlsx (save as B3_METH_CELL_DATA_SampleAnnotations.csv)
• Cancer type / cell line annotation: http://www.cancerrxgene.org/gdsc1000/GDSC1000_WebResources//Data/suppData/TableS1E.xlsx (remove row 1, 2 and 4; remove column 1; and save as sampleDesc.csv)
• Gene expression: http://www.cancerrxgene.org/gdsc1000/GDSC1000_WebResources//Data/preprocessed/Cell_line_RMA_proc_basalExp.txt.zip
• Proteomics: http://tcpaportal.org/mclp/download/MCLP-v1.1-Level4.zip (unzip and rename to prot.tsv)
• IC50s: http://www.cancerrxgene.org/gdsc1000/GDSC1000_WebResources//Data/suppData/TableS4A.xlsx (remove first five rows; remove column 1; and save as ic50.csv)

my.merge = function(x, cell_lines) {
  x = merge(x=cell_lines[,1:2], y=x, by.x="COSMIC.identifier", by.y=0)
  rownames(x) = x$Sample.Name
  x = x[,-(1:2)]
  return(x)
}

read.BEM = function(file, cell_lines, transpose=T, data_type) {
  x = read.delim(file, sep="\t", stringsAsFactors=F, row.names=1)
  if(transpose) {
    x = t(x)
  }
  rownames(x) = gsub("^X", "", rownames(x))
  x = my.merge(x, cell_lines)
  
  ind = cell_lines[,data_type]=="Y"
  cell_lines_for_datatype = cell_lines$Sample.Name[ind]
  all_zeros = setdiff(cell_lines_for_datatype, rownames(x))
  x[all_zeros,] = 0
  
  return(x)
}

averageReplicates = function(x, sampleNameCol, exclude) {
  ind = duplicated(x[,sampleNameCol])
  indices = which(ind)
  for(index in indices) {
    ind = which(x[,sampleNameCol]==x[ind,sampleNameCol])
    x[ind[1],-exclude] = colMeans(x[ind,-exclude])
  }
  x = x[-indices,]
  return(x)
}


convert = function(x) {
  new = toupper(gsub("\\.|\\(|\\)|\\-| ","",x))
  return(new)
}

cell_lines = read.csv("data/sampleDesc.csv", stringsAsFactors=F)
rownames(cell_lines) = cell_lines$Sample.Name

mut = read.BEM("data/PANCAN_SEQ_BEM.txt", cell_lines, data_type="Whole.Exome.Sequencing..WES.")
cna = read.BEM("data/PANCAN_CNA_BEM.rdata.txt", cell_lines, transpose=F, data_type="Copy.Number.Alterations..CNA.")

meth = read_tsv("data/F2_METH_CELL_DATA.txt", col_names=T)
meth = as.data.frame(meth)
rownames(meth) = meth[,1]
meth = meth[,-1]
meth = t(meth)
ann = read.delim("data/B3_METH_CELL_DATA_SampleAnnotations.csv", sep=",", header=T, stringsAsFactors=F)
ann$id = paste0(ann$Sentrix_ID, "_", ann$Sentrix_Position)
meth = merge(ann[,c("id","cosmic_id","Sample_Name")], meth, by.x="id", by.y=0)
meth = merge(x=cell_lines[,1:2], y=meth, by.x="COSMIC.identifier", by.y=2)
meth = averageReplicates(meth, "Sample.Name", 1:4)
rownames(meth) = meth$Sample.Name
meth = meth[,-(1:4)]

cancer_type = cell_lines$Cancer.Type..matching.TCGA.label.
names(cancer_type) = cell_lines$Sample.Name
cancer_type[cancer_type=="" | cancer_type=="UNABLE TO CLASSIFY"] = "OTHER"
cancer_type = model.matrix(~cancer_type+0)
colnames(cancer_type) = gsub("^cancer_type","",colnames(cancer_type))
rownames(cancer_type) = cell_lines$Sample.Name

expr = read_tsv("data/Cell_line_RMA_proc_basalExp-1.txt", col_names=T)
expr = as.data.frame(expr)
colnames(expr) = gsub("DATA\\.", "", colnames(expr))
ind = is.na(expr[,1])
expr = expr[!ind,]
rownames(expr) = expr[,1]
expr = expr[,-(1:2)]
expr = t(expr)
expr = my.merge(expr, cell_lines)

ic50 = read.delim("data/ic50.csv", sep=",", header=T, stringsAsFactors=F, row.names=1)

prot = read.delim("data/prot.tsv", sep="\t", stringsAsFactors=F, header=T, row.names=1)
prot = prot[rownames(prot)!="KMH2",] # because we can't determine whether this matches KMH-2 or KM-H2 in GDSC1000
mapping = merge(x=data.frame(prot=rownames(prot), key=convert(rownames(prot))), 
                y=data.frame(gdsc=cell_lines$Sample.Name, key=convert(cell_lines$Sample.Name)), 
                by.x=2, by.y=2)
prot = prot[as.character(mapping$prot),]
rownames(prot) = mapping$gdsc

data = list()
data$mut = mut
data$cna = cna
data$meth = meth
data$cancer_type = cancer_type
data$expr = expr
data$prot = prot
data$ic50 = ic50
data = lapply(data, as.matrix)

Imputation

For proteomics and IC50 data, the data is imputed as described in the manuscript.

trim.nas = function(x) {
  ind1 = colMeans(is.na(x))<0.3
  ind2 = rowMeans(is.na(x[,ind1]))<0.3
  x = x[ind2,ind1]
  return(x)
}

my.impute = function(x) {
  imputed = impute.svd(x)
  imputed = imputed$x
  rownames(imputed) = rownames(x)
  colnames(imputed) = colnames(x)
  return(imputed)
}

data = intersect.samples(data)
data$prot = trim.nas(data$prot)
data = intersect.samples(data)
data$ic50 = trim.nas(data$ic50)
data = intersect.samples(data)
data$prot = my.impute(data$prot)
data$ic50 = my.impute(data$ic50)

Analysis

similarity_funs = replicate(length(data), inner.product)
names(similarity_funs) = names(data)
similarity_funs$cna = similarity_funs$cancer_type = similarity_funs$mut = jaccard
config_matrices = compute.config.matrices(data, similarity_funs)
cors = rv.cor.matrix(config_matrices)

set.seed(1)
nperm = nboots = 1000
cors_boot = run.bootstraps(config_matrices, nperm)
cors_perm = run.permutations(config_matrices, nboots)

suffStat = list(cors=cors, cors_perm=cors_perm)
pc.fit = pc(suffStat, indepTest=rv.link.significance, labels=names(data), alpha=0.01, conservative=T, maj.rule=F, solve.confl=T)

min.pcors.per.edge = function(cors, pc.fit) {
  n = nrow(cors)
  pcors = matrix(NA, n, n)
  for(i in 1:n) {
    for(j in 1:n) {
      if(j %in% pc.fit@graph@edgeL[[i]]$edges) {
        res = c()
        count = 1
        set = setdiff(1:nrow(cors), c(i, j))
        for(k in 1:length(set)) {
          sets = combn(set,k)
          for(l in 1:ncol(sets)) {
            res[count] = rv.pcor(cors, i, j, sets[,l])
            count = count+1
          }
        }
        pcors[i,j] = min(res)
      } else {
        pcors[i,j] = 0
      }
    }
  }
  return(pcors)
}

plot.result = function(pcors, seed=NULL) {
  if(!is.null(seed)) {
    set.seed(seed)
  }
  adj_matrix = pcors
  net = graph_from_adjacency_matrix(adj_matrix, weighted=T)
  E(net)$width = E(net)$weight*10
  E(net)$arrow.size = 0
  V(net)$shape = "none"
  plot(net, edge.color="black")
}

pcors = min.pcors.per.edge(cors, pc.fit)
rownames(pcors) = colnames(pcors) = c("Mutation", "CNA", "Methylation", "Cancer type", "Gene expression", "Proteomics", "Drug response")
plot.result(pcors, seed=11)

get.statistics = function(cors, cors_perm, cors_boot, a, b, set) {
  rv = rv.pcor(cors, a, b, set)
  conf.int = rv.conf.interval(cors_boot, a, b, set, .99)
  p.value = rv.pval(cors, cors_perm, a, b, set)
  statistics = c(rv, conf.int, p.value)
  names(statistics) = c("RV", "conf.interval.low", "conf.interval.high", "p.value")
  return(statistics)
}

a1 = list()
a1[["RV(mutation, IC50)"]] = get.statistics(cors, cors_perm, cors_boot, "mut","ic50",c())
a1[["RV(CNA, IC50)"]] = get.statistics(cors, cors_perm, cors_boot, "cna","ic50",c())
a1[["RV(methylation, IC50)"]] = get.statistics(cors, cors_perm, cors_boot, "meth","ic50",c())
a1[["RV(cancer type, IC50)"]] = get.statistics(cors, cors_perm, cors_boot, "cancer_type","ic50",c())
a1[["RV(mutation, IC50 | expression)"]] = get.statistics(cors, cors_perm, cors_boot, "mut","ic50","expr")
a1[["RV(CNA, IC50 | expression)"]] = get.statistics(cors, cors_perm, cors_boot, "cna","ic50","expr")
a1[["RV(methylation, IC50 | expression)"]] = get.statistics(cors, cors_perm, cors_boot, "meth","ic50","expr")
a1[["RV(cancer type, IC50 | expression)"]] = get.statistics(cors, cors_perm, cors_boot, "cancer_type","ic50","expr")
a1[["RV(mutation, IC50 | proteomics)"]] = get.statistics(cors, cors_perm, cors_boot, "mut","ic50","prot")
a1[["RV(CNA, IC50 | proteomics)"]] = get.statistics(cors, cors_perm, cors_boot, "cna","ic50","prot")
a1[["RV(methylation, IC50 | proteomics)"]] = get.statistics(cors, cors_perm, cors_boot, "meth","ic50","prot")
a1[["RV(cancer type, IC50 | proteomics)"]] = get.statistics(cors, cors_perm, cors_boot, "cancer_type","ic50","prot")
a1 = do.call(rbind, a1)

a2 = list()
a2[["RV(mutation, proteomics)"]]= get.statistics(cors, cors_perm, cors_boot, "mut","prot",c())
a2[["RV(CNA, proteomics)"]] = get.statistics(cors, cors_perm, cors_boot, "cna","prot",c())
a2[["RV(methylation, proteomics)"]] = get.statistics(cors, cors_perm, cors_boot, "meth","prot",c())
a2[["RV(cancer type, proteomics)"]] = get.statistics(cors, cors_perm, cors_boot, "cancer_type","prot",c())
a2[["RV(mutation, proteomics | expression)"]] = get.statistics(cors, cors_perm, cors_boot, "mut","prot","expr")
a2[["RV(CNA, proteomics | expression)"]] = get.statistics(cors, cors_perm, cors_boot, "cna","prot","expr")
a2[["RV(methylation, proteomics | expression)"]] = get.statistics(cors, cors_perm, cors_boot, "meth","prot","expr")
a2[["RV(cancer type, proteomics | expression)"]] = get.statistics(cors, cors_perm, cors_boot, "cancer_type","prot","expr")
a2 = do.call(rbind, a2)

a3 = list()
a3[["RV(expression, IC50)"]] = get.statistics(cors, cors_perm, cors_boot, "expr","ic50",c())
a3[["RV(proteomics, IC50)"]] = get.statistics(cors, cors_perm, cors_boot, "prot","ic50",c())
a3[["RV(expression, IC50 | proteomics)"]] = get.statistics(cors, cors_perm, cors_boot, "expr","ic50","prot")
a3[["RV(proteomics, IC50 | expression)"]] = get.statistics(cors, cors_perm, cors_boot, "prot","ic50","expr")
a3 = do.call(rbind, a3)

my.barplot = function(a, ylim){
  bp = barplot(a[,1], las=2, ylab="(Partial) matrix correlation", ylim=ylim, col="black")
  arrows(bp, a[,2], bp, a[,3], angle=90, code=3, length=0.05, col="blue")
}

ylim = c(min(a1[,2:3], a2[,2:3], a3[,2:3]), max(a1[,2:3], a2[,2:3], a3[,2:3]))

my.barplot(a1, ylim)

my.barplot(a2, ylim)

my.barplot(a3, ylim)

rv.pcor(cors, "expr","prot",c())

## [1] 0.757952

Supplementary Figure 1

Is actually a cartoon that is not based on actual data. Hence, there is no code for Supplementary Figure 1.

Supplementary Figure 2

Drug response prediction analyses. Note that these are more computationally intensive than the rest of analyses. Change the number of cores used according to your machine.

registerDoMC(cores=30)
colnames(data$expr) = paste0(colnames(data$expr),"_expr")
colnames(data$prot) = paste0(colnames(data$prot),"_prot")

x = as.matrix(cbind(data$mut, data$cna, data$meth, data$cancer_type, data$expr, data$prot))
y = data$ic50

x = scale(x, scale=F, center=T)
y = scale(y)

types = c("mut", "cna", "meth", "cancer_type", "expr", "prot")
data_types = c()
for(type in types) {
  data_types = c(data_types, rep(type, ncol(data[[type]])))
}
prot_or_expr = data_types=="prot" | data_types=="expr"

lambda = "lambda.min"
alpha = 0.5

get.perf = function(x, y, method, lambda, alpha, ind=rep(T, ncol(x)), upstream=rep(T, sum(ind))) {
  foreach(i = 1:ncol(y), .combine="c") %dopar% {
    set.seed(1)
    fit = nested.cv(x[,ind], y[,i], upstream, method, alpha=alpha, lambda_glmnet=lambda, lambda_upstream=lambda, lambda_downstream=lambda)
    cor(y[,i], fit$y_hat)
  }
}

perf = list()
perf[["all"]] = get.perf(x, y, "glmnet", lambda, alpha)
perf[["expr+prot"]] = get.perf(x, y, "glmnet", lambda, alpha, ind=prot_or_expr)
perf[["prot_first"]] = get.perf(x, y, "tandem", lambda, alpha, ind=prot_or_expr, upstream=data_types[prot_or_expr]=="prot")
perf[["expr_first"]] = get.perf(x, y, "tandem", lambda, alpha, ind=prot_or_expr, upstream=data_types[prot_or_expr]=="expr")

plot(perf$all, perf$`expr+prot`, pch=16, xlab="Predictive performance all datasets", ylab="Predictive performance gene expression and proteomics")
abline(0,1)

boxplot(perf$all, perf$`expr+prot`, names=c("All datasets", "Gene expression\nand proteomics"), ylab="Predictive performance")

plot(perf$prot_first, perf$expr_first, pch=16, xlab=TeX("Predictive performance $GEX_{unique}$"), ylab=TeX("Predictive performance $PROT_{unique}$"))
abline(0,1)

boxplot(perf$prot_first, perf$expr_first, names=c(TeX("$GEX_{unique}$"), TeX("$PROT_{unique}$")), ylab="Predictive performance")

get.fit = function(x, y, ind=rep(T, ncol(x)), upstream, lambda, alpha) {
  foreach(i = 1:ncol(y)) %dopar% {
    set.seed(1)
    tandem(x[,ind], y[,i], upstream, alpha=alpha, lambda_upstream=lambda, lambda_downstream=lambda)
  }
}

fits = list()
fits[["prot_first"]] = get.fit(x, y, ind=prot_or_expr, upstream=data_types[prot_or_expr]=="prot", lambda, alpha)
fits[["expr_first"]] = get.fit(x, y, ind=prot_or_expr, upstream=data_types[prot_or_expr]=="expr", lambda, alpha)

rel_con = array(NA, dim=c(ncol(y),2,length(fits)))
dimnames(rel_con)[[2]] = c("Gene expression", "Proteomics")
dimnames(rel_con)[[3]] = names(fits)
titles = c(TeX("$GEX_{unique}$"), TeX("$PROT_{unique}$"))
for(type in names(fits)) {
  for(i in 1:ncol(y)) {
    rel_con[i,,type] = relative.contributions(fits[[type]][[i]], x[,prot_or_expr], data_types[prot_or_expr], lambda_glmnet=lambda)
  }
}
for(i in 1:length(fits)) {
  boxplot(rel_con[,,i], ylab="Relative contribution per dataset", pch=16, main=titles[i])
}

get.vi = function(fit, x) {
  beta = coef(fit)[-1]
  y_hat = x%*%beta
  y_hat_ssq = sum(y_hat^2)
  vi = c()
  for(i in 1:ncol(x)) {
    vi[i] = sum((x[,i]*beta[i])^2)
  }
  return(vi)
}

vi = array(NA, dim=c(sum(prot_or_expr),ncol(y),length(fits)))
for(i in 1:length(fits)) {
  for(j in 1:ncol(y)) {
    vi[,j,i] = get.vi(fits[[i]][[j]], x[,prot_or_expr])
  }
}
rownames(vi) = colnames(x)[prot_or_expr]
colnames(vi) = colnames(y)
dimnames(vi)[[3]] = names(fits)

a = apply(vi[,,"prot_first"], 1, mean)
a = sort(a, decreasing=T)
grep("_expr$",names(a), value=T)[1]

## [1] "ABCB1_expr"

a = apply(vi[,,"expr_first"], 1, mean)
a = sort(a, decreasing=T)
grep("_prot$",names(a), value=T)[1:2]

## [1] "MEK1_pS217S221_prot" "EGFR_pY1173_prot"

a = rowMeans(vi[,,"expr_first"])
phospho = grep("_(.)*_prot$",names(a))
abundance = setdiff(grep("_prot$",names(a)), phospho)
beeswarm(list(`Phosphorylation status`=a[phospho],`Protein abundance`=a[abundance]), pch=16, cex=0.8, ylab=TeX("Feature importance $PROT_{unique}$"))

Loaded packages

For reproducability, we include a list of packages and their version numbers.

sessionInfo()

## R version 3.4.3 (2017-11-30)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.4 LTS
## 
## Matrix products: default
## BLAS: /usr/lib/openblas-base/libblas.so.3
## LAPACK: /usr/lib/libopenblasp-r0.2.18.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] RColorBrewer_1.1-2  beeswarm_0.2.3      latex2exp_0.4.0    
##  [4] TANDEM_1.0.2        glmnet_2.0-13       Matrix_1.2-12      
##  [7] doMC_1.3.5          iterators_1.0.9     foreach_1.4.4      
## [10] Rgraphviz_2.20.0    graph_1.54.0        BiocGenerics_0.22.1
## [13] igraph_1.2.1        NMF_0.21.0          cluster_2.0.6      
## [16] rngtools_1.2.4      pkgmaker_0.22       registry_0.5       
## [19] bcv_1.0.1           readr_1.1.1         pcalg_2.5-0        
## [22] iTOP_1.0.0         
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.16      bdsmatrix_1.3-3   lattice_0.20-35  
##  [4] corpcor_1.6.9     rprojroot_1.3-2   digest_0.6.15    
##  [7] gmp_0.5-13.1      V8_1.5            gridBase_0.4-7   
## [10] R6_2.2.2          plyr_1.8.4        backports_1.1.2  
## [13] stats4_3.4.3      evaluate_0.10.1   ggplot2_2.2.1    
## [16] pillar_1.2.1      rlang_0.2.0       lazyeval_0.2.1   
## [19] curl_3.1          rmarkdown_1.9     stringr_1.3.0    
## [22] munsell_0.4.3     compiler_3.4.3    pkgconfig_2.0.1  
## [25] htmltools_0.3.6   tibble_1.4.2      codetools_0.2-15 
## [28] MASS_7.3-47       RBGL_1.52.0       jsonlite_1.5     
## [31] xtable_1.8-2      dagitty_0.2-2     ggm_2.3          
## [34] gtable_0.2.0      magrittr_1.5      scales_0.5.0     
## [37] stringi_1.1.7     reshape2_1.4.3    doParallel_1.0.11
## [40] robustbase_0.92-8 boot_1.3-20       fastICA_1.2-1    
## [43] tools_3.4.3       DEoptimR_1.0-8    sfsmisc_1.1-2    
## [46] hms_0.4.2         abind_1.4-5       clue_0.3-54      
## [49] colorspace_1.3-2  knitr_1.20