SUBTYPES.DATA = list(
list(name='aml', only.primary=F, is.rna.seq=T, is.mirna.seq=T, display.name='AML'),
list(name='breast', only.primary=T, is.rna.seq=T, is.mirna.seq=T, display.name='BIC'),
list(name='colon', only.primary=T, is.rna.seq=T, is.mirna.seq=T, display.name='COAD'),
list(name='gbm', only.primary=T, is.rna.seq=F, is.mirna.seq=F, display.name='GBM'),
list(name='kidney', only.primary=T, is.rna.seq=T, is.mirna.seq=T, display.name='KIRC'),
list(name='liver', only.primary=T, is.rna.seq=T, is.mirna.seq=T, display.name='LIHC'),
list(name='lung', only.primary=T, is.rna.seq=T, is.mirna.seq=T, display.name='LUSC'),
list(name='melanoma', only.primary=F, is.rna.seq=T, is.mirna.seq=T, display.name='SKCM'),
list(name='ovarian', only.primary=T, is.rna.seq=T, is.mirna.seq=T, display.name='OV'),
list(name='sarcoma', only.primary=T, is.rna.seq=T, is.mirna.seq=T, display.name='SARC'))
MAX.NUM.CLUSTERS = 15
OMIC.SUBSETS = list('multi_omics', 'exp', 'methy', 'mirna')
names(OMIC.SUBSETS) = c('all', '1', '2', '3')
get.clustering.results.dir.path <- function() {
return('RESULTS_DIR_PATH')
}
get.plots.dir.path <- function() {
results.dir.path = get.clustering.results.dir.path()
return(file.path(results.dir.path, 'plots'))
}
get.tables.dir.path <- function() {
results.dir.path = get.clustering.results.dir.path()
return(file.path(results.dir.path, 'tables'))
}
subtype.to.display.name <- function(subtype) {
for (i in 1:length(SUBTYPES.DATA)) {
if (SUBTYPES.DATA[[i]]$name == subtype) {
return(SUBTYPES.DATA[[i]]$display.name)
}
}
}
set.omics.list.attr <- function(subtype.raw.data, subtype.data) {
attr(subtype.raw.data[[1]], 'is.seq') = subtype.data$is.rna.seq
attr(subtype.raw.data[[2]], 'is.seq') = F
attr(subtype.raw.data[[3]], 'is.seq') = subtype.data$is.mirna.seq
return(subtype.raw.data)
}
ALGORITHM.NAMES = c('kmeans', 'spectral', 'lracluster', 'pins', 'snf', 'mkl',
'mcca', 'nmf', 'iCluster')
ALGORITHM.DISPLAY.NAMES = as.list(c('K-means', 'Spectral', 'LRAcluster', 'PINS',
'SNF', 'rMKL-LPP', 'MCCA', 'MultiNMF', 'iClusterBayes'))
names(ALGORITHM.DISPLAY.NAMES) = ALGORITHM.NAMES
print.matrix.latex.format <- function(mat) {
print(do.call(paste, as.list(c(colnames(mat), sep=' & '))))
for (i in 1:nrow(mat)) {
print(do.call(paste, as.list(c(rownames(mat)[i], round(mat[i,], digits=2), sep=' & '))))
}
}
perform.all.analyses <- function(benchmark.ret) {
par(mar=rep(1, 4))
for (i in 1:4) {
cur.func = list(benchmark.omics.time, benchmark.omics.num.clusters,
benchmark.omics.surv, benchmark.omics.clinical)[[i]]
for (omic.subset in names(OMIC.SUBSETS)) {
print(paste('current omic subset ', omic.subset))
benchmark.data = cur.func(benchmark.ret, omic.subset)
displayed.benchmark.data = benchmark.data
colnames(displayed.benchmark.data)[1:ncol(displayed.benchmark.data)] =
sapply(as.list(colnames(displayed.benchmark.data)[1:ncol(displayed.benchmark.data)]),
subtype.to.display.name)
rownames(displayed.benchmark.data) = unlist(ALGORITHM.DISPLAY.NAMES[rownames(displayed.benchmark.data)])
print.matrix.latex.format(displayed.benchmark.data)
table.name = c('runtime', 'num_cluster', 'survival', 'clinical')[i]
write.csv(displayed.benchmark.data, file=file.path(get.tables.dir.path(), paste0(table.name, '_', OMIC.SUBSETS[[omic.subset]], '.csv')))
}
print('------------------------')
}
# plots that include all datasets
for (i in 1:4) {
omic.subset = names(OMIC.SUBSETS)[[i]]
benchmark.surv = benchmark.omics.surv(benchmark.ret, omic.subset)
benchmark.clinical = benchmark.omics.clinical(benchmark.ret, omic.subset)
plot.name = list('multi_omics_surv_clinical.tiff', 'exp_surv_clinical.tiff',
'methy_surv_clinical.tiff', 'mirna_surv_clinical.tiff')[[i]]
create.clinical.survival.plots(benchmark.surv, benchmark.clinical, plot.name)
}
# plots for the mean behaviour
create.mean.clinical.survival.plot(benchmark.ret)
}
get.best.single.omic.mat <- function(single.omic.list1, single.omic.list2) {
best.mat1 = matrix(0, ncol=ncol(single.omic.list1[[1]]),
nrow=nrow(single.omic.list1[[1]]))
best.mat2 = matrix(0, ncol=ncol(single.omic.list1[[1]]),
nrow=nrow(single.omic.list1[[1]]))
for (i in 1:nrow(best.mat1)) {
for (j in 1:ncol(best.mat1)) {
values1 = sapply(single.omic.list1, function(mat) mat[i, j])
values2 = sapply(single.omic.list2, function(mat) mat[i, j])
if (any(is.na(values1))) {
best.mat1[i, j] = NA
best.mat2[i, j] = NA
} else {
best.omic = which.max(values1)
best.mat1[i, j] = values1[best.omic]
best.mat2[i, j] = values2[best.omic]
}
}
}
return(list(best.mat1, best.mat2))
}
create.mean.clinical.survival.plot <- function(benchmark.ret) {
tiff(file.path(get.plots.dir.path(), 'mean_surv_clinical.tiff'), width=4500, height=2250, res=300)
layout(matrix(c(1, 1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5), nrow=2, byrow=T))
single.omic.surv.benchmarks = list()
single.omic.clin.benchmarks = list()
for (i in 2:6) {
omic.subset = c(names(OMIC.SUBSETS), 'best_surv', 'best_clin')[i]
alg.cols = c("#8DD3C7", "#BEBADA", "#FB8072", "#80B1D3", "#FDB462", "#B3DE69", "#FCCDE5", "#BC80BD", "#842121")
pch = c(rep(15, 3), rep(3, 3), rep(19, 3))
if (i <= 4) {
benchmark.surv = benchmark.omics.surv(benchmark.ret, omic.subset)
benchmark.clinical = benchmark.omics.clinical(benchmark.ret, omic.subset)
if (i %in% 2:4) {
single.omic.surv.benchmarks[[i - 1]] = benchmark.surv
single.omic.clin.benchmarks[[i - 1]] = benchmark.clinical
}
} else if (i == 5) {
best.mats = get.best.single.omic.mat(single.omic.surv.benchmarks, single.omic.clin.benchmarks)
benchmark.surv = best.mats[[1]]
benchmark.clinical = best.mats[[2]]
} else {
best.mats = get.best.single.omic.mat(single.omic.clin.benchmarks, single.omic.surv.benchmarks)
benchmark.surv = best.mats[[2]]
benchmark.clinical = best.mats[[1]]
}
surv.sum = rowSums(benchmark.surv, na.rm=T)
clin.sum = rowSums(benchmark.clinical, na.rm=T)
# for mcca, the sum is for an empty vector, and is equal 0, so we remove this
available.indices = surv.sum != 0
surv.sum = surv.sum[available.indices]
clin.sum = clin.sum[available.indices]
alg.cols = alg.cols[available.indices]
pch = pch[available.indices]
if (omic.subset == '1') {
xlab = '-log10(logrank pvalue)'
ylab = '# enriched clinical parameters'
} else {
xlab = ''
ylab = ''
}
if (i %in% c(1, 5, 6)) {
y.min = min(clin.sum) - 1
x.min = min(surv.sum) - 1
} else {
y.min = 0
x.min = 0
}
subplot.name = c('Multi-omics', 'Gene Expression', 'DNA Methylation', 'miRNA Expression',
'Best single-omic (survival)', 'Best single-omic (clinical)')[i]
plot(surv.sum, clin.sum, main=subplot.name, xlab=xlab, ylab=ylab,
xlim=c(x.min, max(surv.sum) + 1), ylim=c(y.min, max(clin.sum + 1)), col=alg.cols, pch=pch, cex.lab=1.8, cex=4, cex.axis=1.5, cex.main=2.5, lwd=4)
}
dev.off()
}
create.clinical.survival.plots <- function(benchmark.surv, benchmark.clinical, plot.name) {
num.subtypes = ncol(benchmark.surv)
tiff(file.path(get.plots.dir.path(), plot.name), width=4500, height=1875, res=300)
alg.cols = c("#8DD3C7", "#BEBADA", "#FB8072", "#80B1D3", "#FDB462", "#B3DE69", "#FCCDE5", "#BC80BD", "#842121")
pch = c(rep(15, 3), rep(3, 3), rep(19, 3))
par(mfrow=c(2, num.subtypes / 2), mar=c(4.1, 4, 3.2, 1))
for (i in 1:num.subtypes) {
subtype = colnames(benchmark.surv)[i]
subtype.surv = benchmark.surv[,subtype]
subtype.clinical = benchmark.clinical[,subtype]
available.indices = !is.na(subtype.surv)
subtype.surv = subtype.surv[available.indices]
subtype.clinical = subtype.clinical[available.indices]
current.cols = alg.cols[available.indices]
current.pch = pch[available.indices]
surv.significance = -log10(0.05)
if (i == 1) {
xlab = '-log10(logrank pvalue)'
ylab = '# enriched clinical parameters'
} else {
xlab = ''
ylab = ''
}
plot(subtype.surv, subtype.clinical, main=subtype.to.display.name(subtype), xlab=xlab, ylab=ylab,
xlim=c(0, max(subtype.surv, surv.significance) + 0.2), ylim=c(0, max(subtype.clinical + 1)), col=current.cols, pch=current.pch, cex.lab=1.4, cex=2.4, cex.axis=1.5, cex.main=1.8, lwd=3)
abline(v=surv.significance, col='red')
}
dev.off()
}
plot.legend <- function() {
alg.cols = c("#8DD3C7", "#BEBADA", "#FB8072", "#80B1D3", "#FDB462", "#B3DE69", "#FCCDE5", "#BC80BD", "#842121")
pch = c(rep(15, 3), rep(3, 3), rep(19, 3))
lwds = c(rep(NA, 3), rep(3, 3), rep(NA, 3))
pt.cexs = 1.8
algorithm.names = ALGORITHM.DISPLAY.NAMES
algorithm.names[length(algorithm.names)] = paste0(algorithm.names[length(algorithm.names)], ' ')
width=4200
if (F) {
alg.cols = alg.cols[-7]
pch = pch[-7]
lwds = lwds[-7]
algorithm.names = algorithm.names[-7]
pt.cexs = pt.cexs[-7]
width=3800
}
tiff(file.path(get.plots.dir.path(), 'legend.tif'), width=width, res=300)
par(font=2, mar=rep(1, 4))
plot(0,xaxt='n',yaxt='n',bty='n',pch='',ylab='',xlab='')
legend(0.58, 1, legend=algorithm.names, col=alg.cols, pch=pch, horiz=T, pt.cex=pt.cexs, pt.lwd=lwds)
dev.off()
}
analyze.benchmark <- function() {
all.clusterings = list()
all.timings = list()
for (i in 1:length(SUBTYPES.DATA)) {
current.subtype.data = SUBTYPES.DATA[[i]]
subtype = current.subtype.data$name
subtype.raw.data = get.raw.data(subtype,
only.primary=current.subtype.data$only.primary)
all.clusterings[[subtype]] = list()
all.timings[[subtype]] = list()
for (algorithm.name in ALGORITHM.NAMES) {
all.clusterings[[subtype]][[algorithm.name]] = list()
all.timings[[subtype]][[algorithm.name]] = list()
for (j in c('all', '1', '2', '3')) {
clustering.path = file.path(get.clustering.results.dir.path(),
paste(subtype, algorithm.name, j, sep='_'))
timing.path = file.path(get.clustering.results.dir.path(),
paste(subtype, algorithm.name, j, 'timing', sep='_'))
load(clustering.path)
load(timing.path)
if (!any(is.na(clustering))) {
names(clustering) = colnames(subtype.raw.data[[1]])
}
all.clusterings[[subtype]][[algorithm.name]][[j]] = clustering
all.timings[[subtype]][[algorithm.name]][[j]] = timing
}
}
}
return(list(all.clusterings=all.clusterings, all.timings=all.timings))
}
check.empirical.surv <- function(old.pvals, new.pvals) {
is.in.conf = matrix(0, ncol=ncol(old.pvals), nrow=nrow(old.pvals))
for (i in 1:nrow(old.pvals)) {
for (j in 1:ncol(old.pvals)) {
old.pval = old.pvals[i, j]
if (is.na(old.pval)) {
is.in.conf[i, j] = NA
next
}
if (old.pval == 1e-10) old.pval = 1e-5
num.runs = floor(30 / old.pval)
new.pval = new.pvals[i, j]
num.success = round(new.pval * num.runs)
print(c(num.success, num.runs))
conf.int = binom.test(num.success, num.runs)$conf.int
in.conf.int = old.pval >= conf.int[1] & old.pval <= conf.int[2]
is.in.conf[i, j] = in.conf.int
}
}
return(is.in.conf)
}
get.empirical.surv <- function(clustering, subtype) {
set.seed(42)
surv.ret = check.survival(clustering, subtype)
orig.chisq = surv.ret$chisq
orig.pvalue = get.logrank.pvalue(surv.ret)
# The initial number of permutations to run
num.perms = round(min(max(10 / orig.pvalue, 1000), 1e6))
should.continue = T
total.num.perms = 0
total.num.extreme.chisq = 0
while (should.continue) {
print('Another iteration in empirical survival calculation')
print(num.perms)
perm.chisq = as.numeric(mclapply(1:num.perms, function(i) {
cur.clustering = sample(clustering)
names(cur.clustering) = names(clustering)
cur.chisq = check.survival(cur.clustering, subtype)$chisq
return(cur.chisq)
}, mc.cores=50))
total.num.perms = total.num.perms + num.perms
total.num.extreme.chisq = total.num.extreme.chisq + sum(perm.chisq >= orig.chisq)
binom.ret = binom.test(total.num.extreme.chisq, total.num.perms)
cur.pvalue = binom.ret$estimate
cur.conf.int = binom.ret$conf.int
print(c(total.num.extreme.chisq, total.num.perms))
print(cur.pvalue)
print(cur.conf.int)
sig.threshold = 0.05
is.conf.small = ((cur.conf.int[2] - cur.pvalue) < min(cur.pvalue / 10, 0.01)) & ((cur.pvalue - cur.conf.int[1]) < min(cur.pvalue / 10, 0.01))
is.threshold.in.conf = cur.conf.int[1] < sig.threshold & cur.conf.int[2] > sig.threshold
if ((is.conf.small & !is.threshold.in.conf) | (total.num.perms > 2e7)) {
#if (is.conf.small) {
should.continue = F
} else {
num.perms = 1e5
}
}
return(list(pvalue = cur.pvalue, conf.int = cur.conf.int, total.num.perms=total.num.perms,
total.num.extreme.chisq=total.num.extreme.chisq))
}
benchmark.omics.surv <- function(benchmark.results, omics='all') {
all.surv.pvalues = matrix(1, ncol=length(SUBTYPES.DATA), nrow=length(ALGORITHM.NAMES))
rownames(all.surv.pvalues) = ALGORITHM.NAMES
colnames(all.surv.pvalues) = sapply(SUBTYPES.DATA, function(x) x$name)
all.clusterings = benchmark.results$all.clusterings
for (i in 1:length(all.clusterings)) {
subtype = colnames(all.surv.pvalues)[i]
subtype.clusterings = all.clusterings[[subtype]]
for (j in 1:length(subtype.clusterings)) {
surv.path = file.path(get.clustering.results.dir.path(),
paste(subtype, ALGORITHM.NAMES[j], omics, 'surv', sep='_'))
if (file.exists(surv.path)) {
load(surv.path)
pvalue = empirical.surv.ret$pvalue
} else {
clustering = subtype.clusterings[[ALGORITHM.NAMES[j]]][[omics]]
if (length(table(clustering)) > 1) {
#pvalue = -log10(get.logrank.pvalue(check.survival(clustering, subtype)))
empirical.surv.ret = get.empirical.surv(clustering, subtype)
save(empirical.surv.ret, file=surv.path)
pvalue = empirical.surv.ret$pvalue
} else {
pvalue = NA
}
}
all.surv.pvalues[j, i] = -log10(pvalue)
}
}
return(all.surv.pvalues)
}
benchmark.omics.num.clusters <- function(benchmark.results, omics='all') {
num.clusters = matrix(1, ncol=length(SUBTYPES.DATA), nrow=length(ALGORITHM.NAMES))
rownames(num.clusters) = ALGORITHM.NAMES
colnames(num.clusters) = sapply(SUBTYPES.DATA, function(x) x$name)
all.clusterings = benchmark.results$all.clusterings
for (i in 1:length(all.clusterings)) {
subtype = colnames(num.clusters)[i]
subtype.clusterings = all.clusterings[[subtype]]
for (j in 1:length(subtype.clusterings)) {
clustering = subtype.clusterings[[ALGORITHM.NAMES[j]]][[omics]]
num.clusters[j, i] = max(clustering)
}
}
return(num.clusters)
}
benchmark.omics.clinical <- function(benchmark.results, omics='all') {
num.clinical.enrich = matrix(1, ncol=length(SUBTYPES.DATA), nrow=length(ALGORITHM.NAMES))
rownames(num.clinical.enrich) = ALGORITHM.NAMES
colnames(num.clinical.enrich) = sapply(SUBTYPES.DATA, function(x) x$name)
total.num.tested.parameters = 0
all.clusterings = benchmark.results$all.clusterings
for (i in 1:length(all.clusterings)) {
subtype = colnames(num.clinical.enrich)[i]
subtype.clusterings = all.clusterings[[subtype]]
for (j in 1:length(subtype.clusterings)) {
print('checking clinical enrichment')
print(c(i, j))
clustering = subtype.clusterings[[ALGORITHM.NAMES[j]]][[omics]]
if (any(is.na(clustering))) {
num.clinical.enrich[j, i] = NA
} else {
clin.path = file.path(get.clustering.results.dir.path(),
paste(subtype, ALGORITHM.NAMES[j], omics, 'clin', sep='_'))
if (file.exists(clin.path)) {
load(clin.path)
} else {
enrichment.pvalues = check.clinical.enrichment(clustering, subtype)
save(enrichment.pvalues, file=clin.path)
}
if (j == 1) {
total.num.tested.parameters = total.num.tested.parameters + length(enrichment.pvalues)
}
num.clinical.enrich[j, i] = sum(enrichment.pvalues * length(enrichment.pvalues) < 0.05)
}
}
}
print(paste0('Total number of parameters tested:', total.num.tested.parameters))
return(num.clinical.enrich)
}
benchmark.omics.time <- function(benchmark.results, omics='all') {
all.alg.times = matrix(1, ncol=length(SUBTYPES.DATA), nrow=length(ALGORITHM.NAMES))
rownames(all.alg.times) = ALGORITHM.NAMES
colnames(all.alg.times) = sapply(SUBTYPES.DATA, function(x) x$name)
all.timings = benchmark.results$all.timings
for (i in 1:length(all.timings)) {
subtype = colnames(all.alg.times)[i]
subtype.timings = all.timings[[subtype]]
for (j in 1:length(subtype.timings)) {
timing = subtype.timings[[ALGORITHM.NAMES[j]]][[omics]]
all.alg.times[j, i] = timing
}
}
return(all.alg.times)
}
benchmark.alg.agreement <- function(benchmark.results, omics='all') {
all.clusterings = benchmark.results$all.clusterings
rand.matrix = matrix(NA, ncol=length(ALGORITHM.NAMES), nrow=length(ALGORITHM.NAMES))
colnames(rand.matrix) = ALGORITHM.NAMES
rownames(rand.matrix) = ALGORITHM.NAMES
for (alg1.index in 1:length(ALGORITHM.NAMES)) {
alg1 = ALGORITHM.NAMES[alg1.index]
for (alg2.index in 1:length(ALGORITHM.NAMES)) {
alg2 = ALGORITHM.NAMES[alg2.index]
rands = c()
for (i in 1:length(SUBTYPES.DATA)) {
clustering1 = all.clusterings[[i]][[alg1]][[omics]]
clustering2 = all.clusterings[[i]][[alg2]][[omics]]
rands = c(rands, adj.rand.index(clustering1, clustering2))
}
mean.rand = mean(rands)
rand.matrix[alg1.index, alg2.index] = mean.rand
}
}
return(rand.matrix)
}
get.logrank.pvalue <- function(survdiff.res) {
1 - pchisq(survdiff.res$chisq, length(survdiff.res$n) - 1)
}
run.benchmark <- function() {
for (i in 1:length(SUBTYPES.DATA)) {
current.subtype.data = SUBTYPES.DATA[[i]]
subtype = current.subtype.data$name
subtype.raw.data = get.raw.data(subtype,
only.primary=current.subtype.data$only.primary)
subtype.raw.data = set.omics.list.attr(subtype.raw.data,
current.subtype.data)
for (algorithm.name in ALGORITHM.NAMES) {
for (j in c('all', '1', '2', '3')) {
set.seed(42)
print(paste('subtype', subtype, 'running algorithm', algorithm.name, j))
clustering.path = file.path(get.clustering.results.dir.path(),
paste(subtype, algorithm.name, j, sep='_'))
timing.path = file.path(get.clustering.results.dir.path(),
paste(subtype, algorithm.name, j, 'timing', sep='_'))
if (!file.exists(clustering.path)) {
algorithm.func.name = paste0('run.', algorithm.name)
algorithm.func = get(algorithm.func.name)
if (j == 'all') {
cur.iteration.data = subtype.raw.data
} else {
cur.iteration.data = subtype.raw.data[as.numeric(j)]
}
algorithm.ret = algorithm.func(cur.iteration.data, current.subtype.data)
clustering = algorithm.ret$clustering
timing = algorithm.ret$timing
print('before saving')
save(clustering, file = clustering.path)
save(timing, file = timing.path)
}
}
}
}
}
get.dataset.dir.path <- function() {
return('DATASETS_PATH')
}
log.and.normalize <- function(omics.data, subtype.data, normalize=T,
filter.var=F) {
# filter features with no variance at all
for (i in 1:length(omics.data)) {
omics.data[[i]] = omics.data[[i]][apply(omics.data[[i]], 1, var) > 0,]
}
for (i in 1:length(omics.data)) {
if (attr(omics.data[[i]], 'is.seq')) {
omics.data[[i]] = log(1+omics.data[[i]])
}
}
if (filter.var) {
omics.data = lapply(omics.data, keep.high.var.features)
}
if (normalize) {
omics.data = lapply(omics.data, normalize.matrix)
}
return(omics.data)
}
normalize.matrix <- function(data.matrix) {
temp = data.matrix - rowMeans(data.matrix)
should.keep = (apply(temp, 1, sd) != 0)
return ((temp / apply(temp, 1, sd))[should.keep, ])
}
filter.non.tumor.samples <- function(raw.datum, only.primary=only.primary) {
# 01 is primary, 06 is metastatic, 03 is blood derived cancer
if (!only.primary)
return(raw.datum[,substring(colnames(raw.datum), 14, 15) %in% c('01', '03', '06')])
else
return(raw.datum[,substring(colnames(raw.datum), 14, 15) %in% c('01')])
}
get.fixed.names <- function(patient.names, include.type=F) {
# fix the TCGA names to only include the patient ids
if (include.type) {
return(gsub('-', '\\.', toupper(substring(patient.names, 1, 15))))
} else {
return(gsub('-', '\\.', toupper(substring(patient.names, 1, 12))))
}
}
fix.patient.names <- function(subtype.raw.data, include.type=F) {
for (i in 1:length(subtype.raw.data)) {
colnames(subtype.raw.data[[i]]) = get.fixed.names(colnames(subtype.raw.data[[i]]),
include.type)
}
return(subtype.raw.data)
}
get.raw.data <- function(subtype.name,
datasets.path = get.dataset.dir.path(),
only.primary=NA) {
omics.dir = file.path(datasets.path, subtype.name)
omics.files = list.files(omics.dir)
omics.files = setdiff(omics.files, c('survival'))
raw.data = lapply(file.path(omics.dir, omics.files), read.table)
if (!is.na(only.primary)) {
raw.data = lapply(raw.data, function(x) filter.non.tumor.samples(x, only.primary = only.primary))
}
name.corrected.data = fix.patient.names(raw.data)
patients.intersection = Reduce(intersect, lapply(name.corrected.data, colnames))
ret.data = lapply(name.corrected.data, function(datum) datum[,patients.intersection])
return(ret.data)
}
get.elbow <- function(values, is.max) {
second.derivatives = c()
for (i in 2:(length(values) - 1)) {
second.derivative = values[i + 1] + values[i - 1] - 2 * values[i]
second.derivatives = c(second.derivatives, second.derivative)
}
print(second.derivatives)
if (is.max) {
return(which.max(second.derivatives) + 1)
} else {
return(which.min(second.derivatives) + 1)
}
}
# Does not support a single omic dataset
run.mcca <- function(omics.list, subtype.data) {
if (length(omics.list) == 1) {
return(list(clustering=rep(NA, ncol(omics.list[[1]])), timing=1))
}
start = Sys.time()
omics.list = log.and.normalize(omics.list, subtype.data,
normalize = T,
filter.var = T)
subtype = subtype.data$name
omics.transposed = lapply(omics.list, t)
cca.ret = PMA::MultiCCA(omics.transposed,
ncomponents = MAX.NUM.CLUSTERS)
sample.rep = omics.transposed[[1]] %*% cca.ret$ws[[1]]
explained.vars = sapply(1:MAX.NUM.CLUSTERS,
function(i) sum(unlist(apply(sample.rep[1:i,,drop=F], 2, var))))
dimension = get.elbow(explained.vars, is.max=F)
print(dimension)
sample.rep = sample.rep[,1:dimension]
sils = c()
clustering.per.num.clusters = list()
for (num.clusters in 2:MAX.NUM.CLUSTERS) {
cur.clustering = kmeans(sample.rep, num.clusters, iter.max=100, nstart=30)$cluster
sil = get.clustering.silhouette(list(t(sample.rep)), cur.clustering)
sils = c(sils, sil)
clustering.per.num.clusters[[num.clusters - 1]] = cur.clustering
}
# NOTE: the next line contains an error. We mistakenly selected the minimal rather maximal silhouette.
# See more details in: http://acgt.cs.tau.ac.il/multi_omic_benchmark/download.html.
cca.clustering = clustering.per.num.clusters[[which.min(sils)]]
time.taken = as.numeric(Sys.time() - start, units='secs')
return(list(clustering=cca.clustering, timing=time.taken))
}
run.snf <- function(omics.list, subtype.data) {
start = Sys.time()
omics.list = log.and.normalize(omics.list, subtype.data)
subtype = subtype.data$name
alpha=0.5
T.val=30
num.neighbors = round(ncol(omics.list[[1]]) / 10)
similarity.data = lapply(omics.list, function(x) {affinityMatrix(dist2(as.matrix(t(x)),as.matrix(t(x))),
num.neighbors, alpha)})
if (length(similarity.data) == 1) {
W = similarity.data[[1]]
} else {
W = SNF(similarity.data, num.neighbors, T.val)
}
num.clusters = estimateNumberOfClustersGivenGraph(W, 2:MAX.NUM.CLUSTERS)[[3]]
clustering = spectralClustering(W, num.clusters)
time.taken = as.numeric(Sys.time() - start, units='secs')
return(list(clustering=clustering, timing=time.taken))
}
run.iCluster <- function(omics.list, subtype.data) {
omics.list = log.and.normalize(omics.list, subtype.data, normalize = F)
start = Sys.time()
subtype = subtype.data$name
dev.ratios = c()
icluster.rets = list()
if (length(omics.list) == 1) {
icluster.ret = iClusterPlus::tune.iClusterBayes(cpus=(MAX.NUM.CLUSTERS - 1), t(omics.list[[1]]),
K=1:(MAX.NUM.CLUSTERS - 1), type=c('gaussian'))$fit
} else {
icluster.ret = iClusterPlus::tune.iClusterBayes(cpus=(MAX.NUM.CLUSTERS - 1), t(omics.list[[1]]),
t(omics.list[[2]]),
t(omics.list[[3]]),
K=1:(MAX.NUM.CLUSTERS - 1), type=rep('gaussian', 3))$fit
}
dev.ratios = lapply(1:(MAX.NUM.CLUSTERS - 1), function(i) icluster.ret[[i]]$dev.ratio)
print('dev.ratios are:')
print(dev.ratios)
optimal.solution = icluster.ret[[which.max(dev.ratios)]]
time.taken = as.numeric(Sys.time() - start, units='secs')
return(list(clustering=optimal.solution$clusters,
timing=time.taken))
}
get.mkl.binary.path = function() {
return('MKL_BINARY_PATH')
}
get.mkl.arguments.path = function() {
return('MKL_ARGS_PATH')
}
run.mkl <- function(omics.list, subtype.data) {
start = Sys.time()
omics.list = log.and.normalize(omics.list, subtype.data)
subtype = subtype.data$name
omics.list = lapply(omics.list, normalize.matrix)
time.taken = as.numeric(Sys.time() - start, units='secs')
export.subtype.to.mkl(omics.list, subtype)
start = Sys.time()
bin.path = get.mkl.binary.path()
subtype.dir = paste0(get.mkl.arguments.path(), subtype, '\\')
paste0(subtype.dir, 'kernels')
command = paste(bin.path, paste0(subtype.dir, 'kernels'),
paste0(subtype.dir, 'ids'),
paste0(subtype.dir, 'output'),
'9', '5')
command.return = system(command)
stopifnot(command.return == 0)
time.taken2 = as.numeric(Sys.time() - start, units='secs')
clustering = get.mkl.clustering(subtype)
return(list(clustering=clustering,
timing=time.taken + time.taken2))
}
run.nmf <- function(omics.list, subtype.data) {
total.time.taken = 0
start = Sys.time()
omics.list = log.and.normalize(omics.list, subtype.data,
filter.var = T, normalize = F)
time.taken = as.numeric(Sys.time() - start, units='secs')
total.time.taken = total.time.taken + time.taken
save.subtype.matlab.format(omics.list)
subtype = subtype.data$name
if (length(omics.list) > 1) {
command.ret = system('MULTI_NMF_COMMAND')
stopifnot(command.ret == 0)
nmf.timing = read.csv('SOME_TEMP_PATH_TIMING', header=F)[1, 1]
total.time.taken = total.time.taken + nmf.timing
} else {
for (k in 1:MAX.NUM.CLUSTERS) {
start = Sys.time()
file.name = paste0('SOME_TEMP_PATH/', k, '_consensus')
nmf.ret = nmf(omics.list[[1]], k, method='lee')
coef.mat = t(coef(nmf.ret))
time.taken = as.numeric(Sys.time() - start, units='secs')
total.time.taken = total.time.taken + time.taken
write.table(coef.mat, file=file.name, quote=F, row.names=F, col.names=F, sep=',')
}
}
explained.vars = c()
clustering.per.num.clusters = list()
for (k in 1:MAX.NUM.CLUSTERS) {
file.name = paste0('SOME_TEMP_PATH/', k, '_consensus')
consensus.mat = read.csv(file.name, header=F)
start = Sys.time()
cur.clustering = apply(consensus.mat, 1, which.max)
explained.var = sum(unlist(apply(consensus.mat, 2, var)))
explained.vars = c(explained.vars, explained.var)
clustering.per.num.clusters[[k]] = cur.clustering
time.taken = as.numeric(Sys.time() - start, units='secs')
total.time.taken = total.time.taken + time.taken
}
dimension = get.elbow(explained.vars, is.max=F)
nmf.clustering = clustering.per.num.clusters[[dimension]]
return(list(clustering=nmf.clustering, timing=total.time.taken))
}
run.pins <- function(omics.list, subtype.data) {
start = Sys.time()
omics.list = log.and.normalize(omics.list, subtype.data, normalize = F)
subtype = subtype.data$name
omics.transposed = lapply(omics.list, t)
if (length(omics.list) == 1) {
pins.ret = PINSPlus::PerturbationClustering(data=omics.transposed[[1]],
kMax = MAX.NUM.CLUSTERS)
clustering = pins.ret$cluster
} else {
pins.ret = PINSPlus::SubtypingOmicsData(dataList=omics.transposed,
kMax = MAX.NUM.CLUSTERS)
clustering = pins.ret$cluster2
}
time.taken = as.numeric(Sys.time() - start, units='secs')
return(list(clustering=clustering, timing=time.taken))
}
run.lracluster <- function(omics.list, subtype.data) {
omics.list = log.and.normalize(omics.list, subtype.data, normalize = F)
subtype = subtype.data$name
start = Sys.time()
dim.range = 1:MAX.NUM.CLUSTERS
all.clustering.results = list()
omics.matrix.list = lapply(omics.list, as.matrix)
for (dimension in dim.range) {
print(paste('running lra cluster for dimension', dimension))
data.names = c('gene expression', 'methylation', 'miRNA expression')
clustering.results = LRAcluster(omics.matrix.list,
rep('gaussian', length(omics.list)),
dimension=dimension, data.names)
all.clustering.results[[dimension]] = clustering.results
}
explained.var = sapply(all.clustering.results, function(x) x$potential)
print(explained.var)
dimension = get.elbow(explained.var, is.max=F)
print(dimension)
solution = all.clustering.results[[dimension]]$coordinate
sils = c()
clustering.per.num.clusters = list()
for (num.clusters in 2:MAX.NUM.CLUSTERS) {
print(paste('running kmeans in lra cluster for num clusters', num.clusters))
cur.clustering = kmeans(t(solution), num.clusters, iter.max=100, nstart=60)$cluster
sil = get.clustering.silhouette(list(solution), cur.clustering)
sils = c(sils, sil)
clustering.per.num.clusters[[num.clusters - 1]] = cur.clustering
}
print(sils)
# NOTE: the next line contains an error. We mistakenly selected the minimal rather maximal silhouette.
# See more details in: http://acgt.cs.tau.ac.il/multi_omic_benchmark/download.html.
chosen.clustering = clustering.per.num.clusters[[which.min(sils)]]
time.taken = as.numeric(Sys.time() - start, units='secs')
return(list(clustering=chosen.clustering, timing=time.taken))
}
run.kmeans <- function(omics.list, subtype.data) {
start = Sys.time()
omics.list = log.and.normalize(omics.list, subtype.data,
filter.var = T)
subtype = subtype.data$name
all.withinss = c()
all.clusterings = list()
k.range = 1:MAX.NUM.CLUSTERS
for (k in k.range) {
concat.omics = do.call(rbind, omics.list)
kmeans.ret = kmeans(t(concat.omics), k, iter.max=100, nstart=60)
all.withinss = c(all.withinss, kmeans.ret$tot.withinss)
all.clusterings[[k]] = kmeans.ret$cluster
}
best.k = get.elbow(all.withinss, is.max=T)
time.taken = as.numeric(Sys.time() - start, units='secs')
return(list(clustering=all.clusterings[[best.k]],
timing=time.taken))
}
run.spectral <- function(omics.list, subtype.data) {
start = Sys.time()
omics.list = log.and.normalize(omics.list, subtype.data,
filter.var = T)
subtype = subtype.data$name
concat.omics = do.call(rbind, omics.list)
similarity.data = affinityMatrix(dist2(as.matrix(t(concat.omics)),
as.matrix(t(concat.omics))),
20, 0.5)
num.clusters = estimateNumberOfClustersGivenGraph(similarity.data,
2:MAX.NUM.CLUSTERS)[[3]]
clustering = spectralClustering(similarity.data, num.clusters)
time.taken = as.numeric(Sys.time() - start, units='secs')
return(list(clustering=clustering, timing=time.taken))
}
load.libraries <- function() {
library('PMA')
library('R.matlab')
library('SNFtool')
library('PINSPlus')
#library('LRAcluster')
library('kernlab')
library('survival')
library('NMF')
# bioconductor packages
source("https://bioconductor.org/biocLite.R")
biocLite("impute")
#biocLite("iClusterPlus")
}
########################################
############### MKL ###############
########################################
radial.basis <- function(mat, gamma) {
if (missing(gamma)) {
gamma = 1 / (2*nrow(mat)**2)
}
npatients = ncol(mat)
output.mat = matrix(0, ncol=npatients, nrow=npatients)
for (i in 1:npatients) {
for (j in 1:npatients) {
output.mat[i, j] = exp(-norm(as.matrix(mat[,i] - mat[,j]), type = 'F')**2 * gamma)
}
}
D = apply(output.mat, 2, sum) / npatients
E = sum(D) / npatients
J = matrix(1, nrow=npatients, ncol=1) %*% D
ret = output.mat - J - t(J) + E * matrix(1, ncol=npatients, nrow=npatients)
ret = diag(1/sqrt(diag(ret))) %*% ret %*% diag(1/sqrt(diag(ret)))
return(ret)
}
clear.dir <- function(dir.path) {
files.in.dir = list.files(dir.path)
for (file.in.dir in files.in.dir) {
full.file.path = file.path(dir.path, file.in.dir)
file.remove(full.file.path)
}
}
export.subtype.to.mkl <- function(omics.list, dir.name) {
folder.path = file.path(get.mkl.arguments.path(), dir.name)
if (!dir.exists(folder.path)) {
dir.create(folder.path)
}
kernels.path = file.path(folder.path, 'kernels')
if (!dir.exists(kernels.path)) {
dir.create(kernels.path)
}
clear.dir(kernels.path)
gammas = 10 ** seq(-6, 6, by=3)
for (i in 1:length(omics.list)) {
for (j in 1:length(gammas)) {
datum = omics.list[[i]]
gamma = gammas[[j]] / (2*nrow(datum)**2)
mat = radial.basis(datum, gamma)
R.matlab::writeMat(file.path(kernels.path, paste(i, '_', j, sep='')), mat=mat)
}
}
output.path = file.path(folder.path, 'output')
if (!dir.exists(output.path)) {
dir.create(output.path)
}
clear.dir(output.path)
write.table(colnames(omics.list[[1]]), file=file.path(folder.path, 'ids'),
quote=F, row.names = F, col.names = F)
}
get.mkl.clustering <- function(dir.name) {
folder.path = file.path(get.mkl.arguments.path(), dir.name)
output.path = file.path(folder.path, 'output')
output.files = list.files(output.path)
clustering = read.csv(file.path(output.path, output.files[grep('clusters', output.files)]))[,2]
return(clustering)
}
check.survival <- function(groups, subtype, survival.file.path) {
if (missing(survival.file.path)) {
survival.file.path = get.subtype.survival.path(subtype)
}
survival.data = read.table(survival.file.path, header = TRUE)
patient.names = names(groups)
patient.names.in.file = as.character(survival.data[, 1])
patient.names.in.file = toupper(gsub('-', '\\.', substring(patient.names.in.file, 1, 12)))
stopifnot(all(patient.names %in% patient.names.in.file))
indices = match(patient.names, patient.names.in.file)
ordered.survival.data = survival.data[indices,]
ordered.survival.data["cluster"] <- groups
ordered.survival.data$Survival[is.na(ordered.survival.data$Survival)] = 0
ordered.survival.data$Death[is.na(ordered.survival.data$Death)] = 0
return(survdiff(Surv(Survival, Death) ~ cluster, data=ordered.survival.data))
}
get.subtype.survival.path <- function(subtype) {
datasets.path = get.dataset.dir.path()
survival.file.path = file.path(datasets.path, subtype, 'survival')
return(survival.file.path)
}
get.clustering.silhouette <- function(raw.data, clustering) {
sils = c()
for (i in 1:length(raw.data)) {
x = raw.data[[i]]
distmatrix = dist2(as.matrix(t(x)),as.matrix(t(x)))
sil = silhouette(clustering, dmatrix = distmatrix)[,3]
sils = c(sils, mean(sil))
}
return(mean(sils))
}
get.clinical.params.dir <- function() {
return('CLINICAL_PARAMS_PATH')
}
get.clinical.params <- function(subtype.name) {
clinical.data.path = paste(get.clinical.params.dir(), subtype.name, sep = '')
clinical.params = read.table(clinical.data.path,
sep='\t', header=T, row.names = 1, stringsAsFactors = F)
rownames.with.duplicates = get.fixed.names(rownames(clinical.params))
clinical.params = clinical.params[!duplicated(rownames.with.duplicates),]
rownames(clinical.params) = rownames.with.duplicates[!duplicated(rownames.with.duplicates)]
return(clinical.params)
}
check.clinical.enrichment <- function(clustering, subtype.name) {
clinical.params = get.clinical.params(subtype.name)
clinical.metadata = list(gender='DISCRETE', age_at_initial_pathologic_diagnosis='NUMERIC',
pathologic_M='DISCRETE', pathologic_N='DISCRETE', pathologic_T='DISCRETE', pathologic_stage='DISCRETE')
pvalues = c()
params.being.tested = c()
for (clinical.param in names(clinical.metadata)) {
if (!(clinical.param %in% colnames(clinical.params))) {
#print(paste0('WARNING: ', clinical.param, ' does not appear for subtype ', subtype.name))
next
}
clinical.values = clinical.params[names(clustering),clinical.param]
is.discrete.param = clinical.metadata[clinical.param] == 'DISCRETE'
is.numeric.param = clinical.metadata[clinical.param] == 'NUMERIC'
stopifnot(is.discrete.param | is.numeric.param)
# skip parameter if many missing values
if (is.numeric.param) {
numeric.entries = !is.na(as.numeric(clinical.values))
if (2 * sum(numeric.entries) < length(clinical.values)) {
#print(paste0('WARNING: skipping on ', clinical.param, ' for subtype ', subtype.name))
next
}
} else {
not.na.entries = !is.na(clinical.values)
should.skip = F
if (2 * sum(not.na.entries) < length(clinical.values)) {
should.skip = T
} else if (length(table(clinical.values[not.na.entries])) == 1) {
should.skip = T
}
if (should.skip) {
#print(paste0('WARNING: skipping on ', clinical.param, ' for subtype ', subtype.name))
next
}
}
params.being.tested = c(params.being.tested, clinical.param)
if (is.discrete.param) {
#clustering.with.clinical = cbind(clustering, clinical.values)
#tbl = table(as.data.frame(clustering.with.clinical[!is.na(clinical.values),]))
#test.res = chisq.test(tbl)
#pvalue = test.res$p.value
pvalue = get.empirical.clinical(clustering[!is.na(clinical.values)], clinical.values[!is.na(clinical.values)], T)
} else if (is.numeric.param) {
#test.res = kruskal.test(as.numeric(clinical.values[numeric.entries]),
# clustering[numeric.entries])
#pvalue = test.res$p.value
pvalue = get.empirical.clinical(clustering[numeric.entries], as.numeric(clinical.values[numeric.entries]), F)
}
pvalues = c(pvalues, pvalue)
}
names(pvalues) = params.being.tested
return(pvalues)
}
get.empirical.clinical <- function(clustering, clinical.values, is.chisq) {
set.seed(42)
if (is.chisq) {
clustering.with.clinical = cbind(clustering, clinical.values)
tbl = table(as.data.frame(clustering.with.clinical))
test.res = chisq.test(tbl)
} else {
test.res = kruskal.test(as.numeric(clinical.values), clustering)
}
orig.pvalue = test.res$p.value
num.iter = 1000
total.num.iters = 0
total.num.extreme = 0
should.continue = T
while (should.continue) {
print('another iteration in empirical clinical')
perm.pvalues = as.numeric(mclapply(1:num.iter, function(i) {
cur.clustering = sample(clustering)
names(cur.clustering) = names(clustering)
if (is.chisq) {
clustering.with.clinical = cbind(cur.clustering, clinical.values)
tbl = table(as.data.frame(clustering.with.clinical))
test.res = chisq.test(tbl)
} else {
test.res = kruskal.test(as.numeric(clinical.values), cur.clustering)
}
cur.pvalue = test.res$p.value
return(cur.pvalue)
}, mc.cores=50))
total.num.iters = total.num.iters + num.iter
total.num.extreme = total.num.extreme + sum(perm.pvalues <= orig.pvalue)
binom.ret = binom.test(total.num.extreme, total.num.iters)
cur.pvalue = binom.ret$estimate
cur.conf.int = binom.ret$conf.int
sig.threshold = 0.05
is.threshold.in.conf = cur.conf.int[1] < sig.threshold & cur.conf.int[2] > sig.threshold
if (!is.threshold.in.conf | total.num.iters > 1e5) {
should.continue = F
}
}
return(cur.pvalue)
}
get.nmf.datasets.dir <- function() {
return('NMF_DATASETS_PATH')
}
save.subtype.matlab.format <- function(subtype.raw.data) {
data.names = c('1', '2', '3')
full.dir.name = get.nmf.datasets.dir()
dir.create(full.dir.name, showWarnings = F)
for (i in 1:length(subtype.raw.data)) {
full.file.name = file.path(full.dir.name, data.names[i])
write.table(subtype.raw.data[[i]], full.file.name)
}
}
keep.high.var.features <- function(omic, num.features=2000) {
if (nrow(omic) < num.features) {
return(omic)
} else {
feature.vars = apply(omic, 1, var)
threshold = feature.vars[order(feature.vars, decreasing = T)][num.features]
return(omic[feature.vars >= threshold,])
}
}
Datasets
Models
