Datasets
Models
Downloads: 1
Diff of
/benchmark.R
[000000]
..
[835948]
Switch to side-by-side view
--- a +++ b/benchmark.R @@ -0,0 +1,1167 @@ +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,]) + } +}