other/nempi_loo.r

In cbg-ethz/nempi: Inferring unobserved perturbations from gene expression data

pipe <- as.character(commandArgs(TRUE)[1])

type <- as.character(commandArgs(TRUE)[2])

bsruns <- as.numeric(commandArgs(TRUE)[3]) # bsruns <- 10

startrun <- as.numeric(commandArgs(TRUE)[4]) # startrun <- 10

if (type %in% c("rf", "svm", "nn")) {
    options(expressions = 500000)
}

safeSave <- function(..., list = character(), file = stop("'file' must be specified"),
                     ascii = FALSE, version = NULL, envir = parent.frame(), compress = isTRUE(!ascii),
                     compression_level, eval.promises = TRUE, precheck = TRUE) {

    i <- 1
    while(file.exists(paste0(file, "_", i))) {
        i <- i + 1
    }

    save(..., list = list, file = paste0(file, "_", i),
         ascii = ascii, version = version, envir = envir, compress = compress,
         compression_level = compression_level, eval.promises = eval.promises, precheck = precheck)
}

aucs <- function(a,b) {
    auc <- roc <- 0
    ppv <- rec <- spec <- NULL
    for (cut in c(2,seq(1,0, length.out = 100),-1)) {
        tmp <- a*0
        tmp[which(a > cut)] <- 1
        tp <- sum(tmp == 1 & b == 1)
        fp <- sum(tmp == 1 & b == 0)
        tn <- sum(tmp == 0 & b == 0)
        fn <- sum(tmp == 0 & b == 1)
        ppvtmp <-  tp/(tp+fp)
        if (is.na(ppvtmp)) { ppvtmp <- 0.5 }
        rectmp <- tp/(tp+fn)
        spectmp <- 1-tn/(tn+fp)
        if (length(ppv) > 0) {
            auc <- auc + (rectmp-rec[length(rec)])*(ppvtmp+ppv[length(ppv)])/2
            roc <- roc + (spectmp-spec[length(spec)])*(rectmp+rec[length(rec)])/2
        }
        ppv <- c(ppv, ppvtmp)
        rec <- c(rec, rectmp)
        spec <- c(spec, spectmp)
    }
    return(list(auc=auc,roc=roc,ppv=ppv,rec=rec,spec=spec))
}

path <- "/cluster/work/bewi/members/mpirkl/"

if (pipe %in% "loo") { 
    library(naturalsort)
    library(nem)
    library(cluster)
    library(Rcpp)
    library(Rgraphviz)
    library(randomForest)
    library(e1071)
    source("mnems.r")
    source("mnems_low.r")
    sourceCpp(code=readChar("mm.cpp", file.info("mm.cpp")$size))
    source("nempi_main.r")
    source("nempi_low.r")
    load("TCGA-BRCA_nempi.rda")
    ## sourceCpp("~/Documents/mnem/src/mm.cpp"); source("~/Documents/mnem/R/mnems.r"); source("~/Documents/mnem/R/mnems_low.r"); source("~/Documents/nempi/R/nempi_main.r"); source("~/Documents/nempi/R/nempi_low.r")
    ## load("nempi_backup/backup3/TCGA-BRCA_nempi.rda")
    P <- Pmut
    P <- P[which(rownames(P) %in% rownames(Rho)), which(colnames(P) %in% colnames(Rho))]
    P <- P[order(rownames(P)), order(colnames(P))]
    P[which(P > 1)] <- 1
    Ptmp <- cbind(P, matrix(0, nrow(P), sum(!(colnames(Rho) %in% colnames(P)))))
    colnames(Ptmp) <- c(colnames(P), colnames(Rho)[which(!(colnames(Rho) %in% colnames(P)))])
    P <- Ptmp[, colnames(Rho)]
    D4 <- D2
    colnames(D4) <- apply(Rho, 2, function(x) {
        Sgenes <- paste(sort(rownames(Rho)[which(x > 0)]), collapse = "_")
        return(Sgenes)
    })
    D5 <- D4[, which(apply(P, 2, sum) != 0)]
    genes.med <- apply(D5, 1, median)
    D5 <- D5[which(genes.med != 0), ]
    PM <- P[, which(apply(P, 2, sum) != 0)]
    tp <- fp <- fn <- tn <- prauc <- numeric(ncol(D5))
    if (type %in% "nempi") {
        converged <- 1
        for (i in startrun:(startrun+bsruns-1)) {
            if (i > ncol(D5)) { break() }
            D6 <- D5[, -i]
            Rho2 <- PM[, -i]
            Rho3 <- PM
            Rho3[, i] <- 0
            nempitmp <- nempi(D6, Gamma = Rho2, full = TRUE, converged = converged)
            Rho4 <- Rho3
            Rho4[, -i] <- nempitmp$Gamma
            nempitmp <- nempi(D5, Gamma = Rho4, full = FALSE, converged = converged)
            Rho5 <- t(mytc(nempitmp$res$adj))%*%nempitmp$Gamma
            tp[i] <- sum(PM[, i] == 1 & Rho5[, i] >= 1/8)
            fp[i] <- sum(PM[, i] == 0 & Rho5[, i] >= 1/8)
            fn[i] <- sum(PM[, i] == 1 & Rho5[, i] < 1/8)
            tn[i] <- sum(PM[, i] == 0 & Rho5[, i] < 1/8)
            prauc[i] <- aucs(Rho5[, i],PM[, i])$auc
        }
        safeSave(prauc, tp, fp, fn, tn, file = paste0(path, "nempi_bs/nempi_loo.rda"))
    } else if (type %in% "svm") {
        for (i in startrun:(startrun+bsruns-1)) {
            if (i > ncol(D5)) { break() }
            D6 <- D5[, -i]
            Rho2 <- PM[, -i]
            Rho3 <- PM
            Rho3[, i] <- 0
            svmres <- classpi(D6, full = TRUE, method = "svm")
            Rho4 <- Rho3
            Rho4[, -i] <- svmres$Gamma
            D7 <- D5
            colnames(D7)[i] <- ""
            svmres <- classpi(D7, full = FALSE, method = "svm")
            Rho5 <- svmres$Gamma
            tp[i] <- sum(PM[, i] == 1 & Rho5[, i] >= 1/8)
            fp[i] <- sum(PM[, i] == 0 & Rho5[, i] >= 1/8)
            fn[i] <- sum(PM[, i] == 1 & Rho5[, i] < 1/8)
            tn[i] <- sum(PM[, i] == 0 & Rho5[, i] < 1/8)
            prauc[i] <- aucs(Rho5[, i],PM[, i])$auc
        }
        safeSave(prauc, tp, fp, fn, tn, file = paste0(path, "nempi_bs/svm_loo.rda"))
    } else if (type %in% "rf") {
        for (i in startrun:(startrun+bsruns-1)) {
            if (i > ncol(D5)) { break() }
            D6 <- D5[, -i]
            Rho2 <- PM[, -i]
            Rho3 <- PM
            Rho3[, i] <- 0
            svmres <- classpi(D6, full = TRUE, method = "randomForest")
            Rho4 <- Rho3
            Rho4[, -i] <- svmres$Gamma
            D7 <- D5
            colnames(D7)[i] <- ""
            svmres <- classpi(D7, full = FALSE, method = "randomForest")
            Rho5 <- svmres$Gamma
            tp[i] <- sum(PM[, i] == 1 & Rho5[, i] >= 1/8)
            fp[i] <- sum(PM[, i] == 0 & Rho5[, i] >= 1/8)
            fn[i] <- sum(PM[, i] == 1 & Rho5[, i] < 1/8)
            tn[i] <- sum(PM[, i] == 0 & Rho5[, i] < 1/8)
            prauc[i] <- aucs(Rho5[, i],PM[, i])$auc
        }
        safeSave(prauc, tp, fp, fn, tn, file = paste0(path, "nempi_bs/rf_loo.rda"))
    } else if (type %in% "nn") {
        for (i in startrun:(startrun+bsruns-1)) {
            if (i > ncol(D5)) { break() }
            D6 <- D5[, -i]
            Rho2 <- PM[, -i]
            Rho3 <- PM
            Rho3[, i] <- 0
            library(nnet)
            svmres <- classpi(D6, full = TRUE, method = "nnet", MaxNWts = 50000, size = 1)
            Rho4 <- Rho3
            Rho4[, -i] <- svmres$Gamma
            D7 <- D5
            colnames(D7)[i] <- ""
            svmres <- classpi(D7, full = FALSE, method = "nnet", MaxNWts = 50000, size = 1)
            Rho5 <- svmres$Gamma
            tp[i] <- sum(PM[, i] == 1 & Rho5[, i] >= 1/8)
            fp[i] <- sum(PM[, i] == 0 & Rho5[, i] >= 1/8)
            fn[i] <- sum(PM[, i] == 1 & Rho5[, i] < 1/8)
            tn[i] <- sum(PM[, i] == 0 & Rho5[, i] < 1/8)
            prauc[i] <- aucs(Rho5[, i],PM[, i])$auc
        }
        safeSave(prauc, tp, fp, fn, tn, file = paste0(path, "nempi_bs/nn_loo.rda"))
    } else if (type %in% "mf") {
        for (i in startrun:(startrun+bsruns-1)) {
            if (i > ncol(D5)) { break() }
            D6 <- D5
            colnames(D6)[i] <- ""
            Rho2 <- PM
            Rho2[, i] <- 0
            Rho3 <- PM
            Rho3[, i] <- 0
            mfdata <- cbind(as.data.frame(t(D6)), factor(colnames(D6)))
            mfdata[which(mfdata == "", arr.ind = TRUE)] <- NA
            library(missForest)
            mfimp <- missForest(mfdata)
            colnames(D6) <- mfimp$ximp[, ncol(mfimp$ximp)]
            tmp <- nem(D6)
            ures <- list()
            ures$Gamma <- apply(getGamma(D6), 2, function(x) return(x/sum(x)))
            ures$res <- list()
            ures$res$adj <- tmp$adj
            ures$null <- TRUE
            ures$combi <- 1
            Rho5 <- ures$Gamma
            tp[i] <- sum(PM[, i] == 1 & Rho5[, i] >= 1/8)
            fp[i] <- sum(PM[, i] == 0 & Rho5[, i] >= 1/8)
            fn[i] <- sum(PM[, i] == 1 & Rho5[, i] < 1/8)
            tn[i] <- sum(PM[, i] == 0 & Rho5[, i] < 1/8)
            prauc[i] <- aucs(Rho5[, i],PM[, i])$auc
        }
        safeSave(prauc, tp, fp, fn, tn, file = paste0(path, "nempi_bs/mf_loo.rda"))
    } else if (type %in% "knn") {
        for (i in startrun:(startrun+bsruns-1)) {
            if (i > ncol(D5)) { break() }
            D6 <- D5
            colnames(D6)[i] <- ""
            Rho2 <- PM
            Rho2[, i] <- 0
            Rho3 <- PM
            Rho3[, i] <- 0
            library(class)
            train <- t(D6[, which(colnames(D6) != "")])
            test <- t(D6[, which(colnames(D6) == "")])
            knn0 <- 0
            cl <- colnames(D6)[which(colnames(D6) != "")]
            knnres <- knn(train, test, cl, prob = TRUE)
            D3 <- D6
            colnames(D3)[which(colnames(D3) %in% "")] <- as.character(knnres)
            tmp <- nem(D3)
            ures <- list()
            ures$Gamma <- apply(getGamma(D3), 2, function(x) return(x/sum(x)))
            ures$res <- list()
            ures$res$adj <- tmp$adj
            ures$null <- TRUE
            ures$combi <- 1
            Rho5 <- ures$Gamma
            tp[i] <- sum(PM[, i] == 1 & Rho5[, i] >= 1/8)
            fp[i] <- sum(PM[, i] == 0 & Rho5[, i] >= 1/8)
            fn[i] <- sum(PM[, i] == 1 & Rho5[, i] < 1/8)
            tn[i] <- sum(PM[, i] == 0 & Rho5[, i] < 1/8)
            prauc[i] <- aucs(Rho5[, i],PM[, i])$auc
        }
        safeSave(prauc, tp, fp, fn, tn, file = paste0(path, "nempi_bs/knn_loo.rda"))
    } else if (type %in% "rand") {
        for (i in startrun:(startrun+bsruns-1)) {
            if (i > ncol(D5)) { break() }
            Rho5 <- PM
            Rho5[,i] <- runif(nrow(Rho5))
            tp[i] <- sum(PM[, i] == 1 & Rho5[, i] >= 1/8)
            fp[i] <- sum(PM[, i] == 0 & Rho5[, i] >= 1/8)
            fn[i] <- sum(PM[, i] == 1 & Rho5[, i] < 1/8)
            tn[i] <- sum(PM[, i] == 0 & Rho5[, i] < 1/8)
            prauc[i] <- aucs(Rho5[, i],PM[, i])$auc
        }
        safeSave(prauc, tp, fp, fn, tn, file = paste0(path, "nempi_bs/rand_loo.rda"))
    }
    stop("successfull loo")
} else {

    library(snowfall)
    library(TCGAbiolinks)

    if (is.null(types)) {
        types <- TCGAbiolinks:::getGDCprojects()$project_id
        ## special "FM-AD"
        dont <- paste0(unique(gsub("-.*", "", types)), "-")
        dont <- dont[-grep("TCGA", dont)]
        samplenr <- matrix(NA, length(types), 2)
        rownames(samplenr) <- types
        donr <- TRUE
        snrcount <- 0
    } else {
        dont <- "AGCT"
        donr <- FALSE
    }

    sizemat <- matrix(0, 1, 2)
    colnames(sizemat) <- c("Tumor", "Normal")
    rownames(sizemat) <- ""

    path <- "mutclust/"

    for (type in types) {
        if (donr) {
            snrcount <- snrcount + 1
        }
        if (length(grep(paste(dont, collapse = "|"), type)) > 0) { next() }
        print(type)
        if (file.exists(paste0(path, type, "_final.rda")) & !newmut & !newllr & !newsave) {
            load(paste0(path, type, "_final.rda"))
        } else {

            summ <- TCGAbiolinks:::getProjectSummary(type)
            library(SummarizedExperiment)

            ## get methylation:
            if (file.exists(paste0(path, type, "_meth.rda"))) {
                load(paste0(path, type, "_meth.rda"))
            } else {
                data <- GDCquery(project = paste(type, sep = ""),
                                 sample.type = "Primary solid Tumor",
                                 data.category = "DNA methylation",
                                 data.type = "Methylation beta value",
                                 legacy = TRUE
                                 )
                GDCdownload(data)
                data <- GDCprepare(data, summarizedExperiment = 0)
                ## map methy sites to genes:
                library(methyAnalysis)
                if (is.data.frame(data)) {
                    meth2 <- as.matrix(data[, -(1:3)])
                    rownames(meth2) <- gsub(";.*", "", data[, 1])
                } else {
                    meth <- data@rowRanges
                    library(TxDb.Hsapiens.UCSC.hg19.knownGene)
                    probe2gene <- annotateDMRInfo(meth, 'TxDb.Hsapiens.UCSC.hg19.knownGene')
                    meth2 <- assay(data)
                    rownames(meth2) <- probe2gene$sigDMRInfo@elementMetadata@listData$GeneSymbol
                }
                meth <- meth2
                meth <- meth[which(apply(meth, 1, function(x) return(any(is.na(x)))) == FALSE), ]
                methm <- meth[which(duplicated(rownames(meth)) == FALSE), ]
                count <- 0
                for (i in which(duplicated(rownames(meth)) == FALSE)) {
                    count <- count + 1
                    methm[count, ] <- apply(meth[which(rownames(meth) %in% rownames(methm)[i]), , drop = FALSE], 2, median)
                }
                meth <- methm
                meth <- meth[order(rownames(meth)), order(colnames(meth))]
                colnames(meth) <- gsub("A$", "", lapply(strsplit(colnames(meth), "-"), function(x) { y <- x[1:4]; y <- paste(y, collapse = "-"); return(y) }))
                methm <- meth[, which(duplicated(colnames(meth)) == FALSE)]
                for (i in which(duplicated(colnames(meth)) == TRUE)) {
                    j <- which(colnames(methm) == colnames(meth)[i])
                    methm[, j] <- apply(meth[, which(colnames(meth) %in% colnames(methm)[i]), drop = FALSE], 2, median)
                }
                meth <- methm
                meth <- meth[order(rownames(meth)), order(colnames(meth))]
                save(meth, file = paste0(path, type, "_meth.rda"))
            }

            print("meth done")

            ## get copy number variation:
            if (file.exists(paste0(path, type, "_cnv.rda"))) {
                load(paste0(path, type, "_cnv.rda"))
            } else {
                data <- getGistic(gsub("TCGA-", "", type), type = "thresholded")
                ## data <- GDCquery(project = paste(type, sep = ""),
                ##                  sample.type = "Primary solid Tumor",
                ##                  data.category = "Copy Number Variation",
                ##                  data.type = "Copy Number Segment",
                ##                  )
                ## GDCdownload(data)
                ## data <- GDCprepare(data)
                cnv <- data[, -(1:3)]
                cnv <- apply(cnv, c(1,2), as.numeric)
                rownames(cnv) <- data[, 1]
                colnames(cnv) <- gsub("A$", "", lapply(strsplit(colnames(cnv), "-"), function(x) { y <- x[1:4]; y <- paste(y, collapse = "-"); return(y) }))
                cnv <- cnv[order(rownames(cnv)), order(colnames(cnv))]
                save(cnv, file = paste0(path, type, "_cnv.rda"))
            }

            print("cnv done")

            ## get expression
            if (file.exists(paste0(path, type, "_query.rda"))) {
                load(paste0(path, type, "_query.rda"))
            } else {
                if (length(grep("-AML$|-LAML$", type)) > 0) {
                    sampletype <- c("Primary Blood Derived Cancer - Peripheral Blood", "Solid Tissue Normal")
                } else {
                    sampletype <- c("Primary Tumor", "Solid Tissue Normal")
                }
                data <- GDCquery(project = paste(type, sep = ""),
                                 sample.type = sampletype,
                                 data.category = "Transcriptome Profiling",
                                 data.type = "Gene Expression Quantification",
                                 workflow.type = "HTSeq - Counts"
                                 )
                GDCdownload(data)
                if (is.null(data)) {
                    print("is null")
                    next()
                }
                data <- GDCprepare(data)
                save(data, file = paste0(path, type, "_query.rda"))
            }
            ## process expression data:
            D <- assay(data)
            class <- data@colData@listData$definition

            print("gene expression done")

            ## get mutation
            if (file.exists(paste0(path, type, "_mut.rda"))) {
                load(paste0(path, type, "_mut.rda"))
            } else {
                type2 <- gsub(paste(paste0(unique(gsub("-.*", "", types)), "-"), collapse = "|"), "", type)
                library(data.table)
                GDCquery_Maf(tumor = type2, save.csv = TRUE, pipeline = "varscan2") # mutation file
                GDCquery_Maf(tumor = type2, save.csv = TRUE, pipeline = "muse") # mutation file
                GDCquery_Maf(tumor = type2, save.csv = TRUE, pipeline = "somaticsniper") # mutation file
                GDCquery_Maf(tumor = type2, save.csv = TRUE, pipeline = "mutect2") # mutation file
                mut <- list()
                count <- 1
                type3 <- gsub("-", "\\.", type)
                for (i in list.files("GDCdata")) {
                    if (length(grep(type3, i)) > 0 & length(grep("csv", i)) > 0) {
                        mut[[count]] <- fread(paste0("GDCdata/", i))
                        count <- count + 1
                    }
                }
                save(mut, file = paste0(path, type, "_mut.rda"))
            }

            ## clinical
            if (file.exists(paste0(path, type, "_clin.rda"))) {
                load(paste0(path, type, "_clin.rda"))
            } else {
                clinical <- GDCquery_clinic(project = type, type = "clinical")
                save(clinical, file = paste0(path, type, "_clin.rda"))
            }

            ## process mutation data: (https://cancer.sanger.ac.uk/cosmic/help/mutation/overview)
            if (file.exists(paste0(path, type, "_mut0.rda")) & !newmut) {
                load(paste0(path, type, "_mut0.rda"))
            } else {
                n <- nmut # try something to get all patients with at least one mutation
                library(snowfall)
                countSamples <- function(i, mut.org, genes, mut.mat, coln) {
                    i <- which(rownames(mut.mat) %in% genes[i])
                    tmp <- mut.org[which(mut.org$Hugo_Symbol %in% rownames(mut.mat)[i]), coln]
                    tmp2 <- mut.mat[i, ]
                    tmp2[which(colnames(mut.mat) %in% tmp)] <- 1
                    return(tmp2)
                }
                typeSamples <- function(i, mut.org, genes, type.mat, coln, coln2) {
                    i <- which(rownames(type.mat) %in% genes[i])
                    tmp <- mut.org[which(mut.org$Hugo_Symbol %in% rownames(type.mat)[i]), coln]
                    tmp3 <- mut.org[which(mut.org$Hugo_Symbol %in% rownames(type.mat)[i]), coln2]
                    tmp2 <- type.mat[i, ]
                    tmp2[which(colnames(type.mat) %in% tmp)] <- tmp3
                    return(tmp2)
                }
                biggl <- list()
                for (i in length(mut)) {
                    mutation <- mut[[i]]
                    biggl[[i]] <- mutation$Hugo_Symbol
                }
                freq <- sort(table(unlist(biggl)), decreasing = TRUE)
                if (n == 0) {
                    allsub <- names(freq)
                } else {
                    allsub <- names(freq)[1:n]
                }
                M <- Mtype <- list()
                for (i in 1:length(mut)) {
                    mutation <- mut[[i]]
                    mut.mat <- matrix(0, length(allsub), length(unique(mutation$Tumor_Sample_Barcode)))
                    type.mat <- matrix("", length(allsub), length(unique(mutation$Tumor_Sample_Barcode)))
                    colnames(type.mat) <- colnames(mut.mat) <- sort(unique(mutation$Tumor_Sample_Barcode))
                    rownames(type.mat) <- rownames(mut.mat) <- allsub
                    coln <- which(colnames(mutation) %in% "Tumor_Sample_Barcode")
                    coln2 <- which(colnames(mutation) %in% "Variant_Classification")
                    mut.org <- mutation[which(mutation$Hugo_Symbol %in% allsub), ]
                    sfInit(parallel = T, cpus = 4)
                    sfExport("mut.mat", "coln")
                    tmp <- sfLapply(as.list(1:length(allsub)), countSamples, mut.org, allsub, mut.mat, coln)
                    sfStop()
                    tmp <- do.call("rbind", tmp)
                    rownames(tmp) <- allsub
                    colnames(tmp) <- colnames(mut.mat)
                    M[[i]] <- tmp
                    sfInit(parallel = T, cpus = 4)
                    sfExport("type.mat", "coln2")
                    tmp <- sfLapply(as.list(1:length(allsub)), typeSamples, mut.org, allsub, type.mat, coln, coln2)
                    sfStop()
                    tmp <- do.call("rbind", tmp)
                    rownames(tmp) <- allsub
                    colnames(tmp) <- colnames(mut.mat)
                    Mtype[[i]] <- tmp
                }
                samples <- intersect(intersect(intersect(colnames(M[[1]]), colnames(M[[2]])), colnames(M[[3]])), colnames(M[[4]]))
                M0 <- M[[1]][, which(colnames(M[[1]]) %in% samples)]
                Mtype0 <- Mtype[[1]][, which(colnames(Mtype[[1]]) %in% samples)]
                for (i in 2:length(M)) {
                    M0 <- M0 + M[[i]][, which(colnames(M[[i]]) %in% samples)]
                    Mtype0 <- matrix(paste(Mtype0, Mtype[[i]][, which(colnames(Mtype[[i]]) %in% samples)]), nrow(Mtype0))
                }
                rownames(Mtype0) <- rownames(M0)
                colnames(Mtype0) <- colnames(M0)
                save(M0, Mtype0, file = paste0(path, type, "_mut0.rda"))
            }

            ## process expression data:
            D <- assay(data)
            class <- data@colData@listData$definition
            M <- M0
            Mtype <- Mtype0
            colnames(M) <- lapply(colnames(M), function(x) {
                y <- unlist(strsplit(x, "-"))
                y <- paste(y[1:4], collapse = "-")
                y <- unlist(strsplit(y, ""))
                y <- paste(y[1:(length(y)-1)], collapse = "")
                return(y)
            })
            colnames(D) <- lapply(colnames(D), function(x) {
                y <- unlist(strsplit(x, "-"))
                y <- paste(y[1:4], collapse = "-")
                y <- unlist(strsplit(y, ""))
                y <- paste(y[1:(length(y)-1)], collapse = "")
                return(y)
            })
            colnames(M) <- gsub("A$", "", lapply(strsplit(colnames(M), "-"), function(x) { y <- x[1:4]; y <- paste(y, collapse = "-"); return(y) }))
            M <- M[order(rownames(M)), order(colnames(M))]
            Mtype <- Mtype[order(rownames(Mtype)), order(colnames(Mtype))]

            print("mutation done")

            ## log odds:
            if (file.exists(paste0(path, type, "_llr.rda")) & !newllr) {
                load(paste0(path, type, "_llr.rda"))
            } else {
                library(edgeR)
                if (sum(class %in% "Solid Tissue Normal") < 10) {
                    distrParNC <- function(i, data) {
                        data[i, ] <- data[i, ] - median(data[, i])
                        llrcol <- numeric(ncol(data))
                        div0 <- quantile(data[i, ], 0.25)
                        div1 <- quantile(data[i, ], 0.75)
                        sigma <- sd(data[i, ])
                        for (j in 1:ncol(data)) {
                            if (data[i, j] <= 0) {
                                llrcol[j] <- log2(div0/data[i, j])
                            } else {
                                llrcol[j] <- log2(data[i, j]/div1)
                            }
                        }
                        return(llrcol)
                    }
                    highcounts <- which(apply(D, 1, median) >= 10)
                    DC <- D[highcounts, ]
                    genenames <- rownames(D)[highcounts, ]
                    nf <- calcNormFactors(DC)
                    DC <- t(t(DC)/nf) # DC <- DC2
                    ## this is very adventurous:
                    ## DC <- t(scale(t(DC)))
                    ## DC <- abs(DC)
                    ## pc <- 1-min(DC)
                    ## tmp <- log2(DC+pc)
                    ## hist(tmp)
                    sfInit(parallel = T, cpus = 4)
                    sfExport("DC")
                    tmp <- do.call("rbind", sfLapply(1:nrow(DC),
                                                     distrParNC, DC))
                    colnames(tmp) <- colnames(DC)
                    rownames(tmp) <- genenames
                    sfStop()
                    DF <- DC
                } else {
                    library(ks)
                    mykcde <- function(x, k) {
                        a <- which.min(abs(k$eval.points - x))
                        b <- k$estimate[a]
                        b <- min(b, 1-b)
                        return(b)
                    }
                    distrParKs <- function(i, data, C) {
                        llrcol <- numeric(ncol(data))
                        ddistr <- list()
                        dogenes <- unique(colnames(data))[which(!(unique(colnames(data)) %in% ""))]
                        for (j in dogenes) {
                            D <- which(colnames(data) %in% j)
                            ddistr[[j]] <- kcde(data[i, D])
                        }
                        cdistr <- kcde(data[i, C])
                        for (j in which(!(colnames(data) %in% ""))) {
                            gene <- colnames(data)[j]
                            llrcol[j] <- log2(mykcde(data[i, j], ddistr[[gene]])/mykcde(data[i,j], cdistr))
                        }
                        llrcol <- llrcol[-C]
                        return(llrcol)
                    }
                    DN <- D[, which(class %in% "Solid Tissue Normal")]
                    DT <- D[, which(class %in% "Primary solid Tumor")]
                    DF <- cbind(DT, DN)
                    C <- (ncol(DT)+1):ncol(DF)
                    highcounts <- which(apply(DF, 1, median) >= 10)
                    DF <- DF[highcounts, ]
                    genenames <- rownames(DF)
                    colnames(DF)[1:ncol(DT)] <- "P" # not knock-down specific!
                    colnames(DF)[grep(paste(unique(gsub("-.*", "", types)), collapse = "|"), colnames(DF))] <- ""
                    nf <- calcNormFactors(DF)
                    DF <- t(t(DF)/nf)
                    sfInit(parallel = T, cpus = 4)
                    sfExport("DF", "C", "mykcde")
                    sfLibrary(ks)
                    tmp <- do.call("rbind", sfLapply(1:nrow(DF),
                                                     distrParKs, DF, C))
                    sfStop()
                    colnames(tmp) <- colnames(DT)
                }
                save(tmp, DF, file = paste0(path, type, "_llr.rda"))
            }

            rownames(tmp) <- rownames(DF)

            par(mfrow=c(1,3))
            hist(tmp)

            tmp[which(is.na(tmp) | is.infinite(tmp))] <- 0

            hist(tmp)

            D <- tmp

            print("expression done")

            ## sd.glob <- sd(tmp)

            ## tmp <- tmp[-which(apply(tmp, 1, sd) < sd.glob), ]

            ## hist(tmp)

            ## prep clinical data:
            clinical[which(clinical$vital_status%in% "dead"), which(colnames(clinical) %in% "vital_status")] <- 1
            clinical[which(clinical$vital_status%in% "alive"), which(colnames(clinical) %in% "vital_status")] <- 0
            count <- 0
            for (stage in sort(unique(clinical$tumor_stage))) {
                clinical[which(clinical$tumor_stage%in% stage), which(colnames(clinical) %in% "tumor_stage")] <- count
                count <- count + 1
            }
            clinical$tumor_stage <- as.numeric(clinical$tumor_stage)
            clinical[which(is.na(clinical$days_to_death)), which(colnames(clinical) %in% "days_to_death")] <- clinical[which(is.na(clinical$days_to_death)), which(colnames(clinical) %in% "days_to_last_follow_up")]
            clinical$vital_status<- as.numeric(clinical$vital_status)

            print("clinical done")

            save(clinical, D, M, Mtype, DF, class, meth, cnv, file = paste0(path, type, "_final.rda"))
        }
        print(table(class))
        sizemat <- rbind(sizemat, table(class))
        rownames(sizemat)[nrow(sizemat)] <- type
        if (donr) {
            samplenr[snrcount, 1] <- sum(class %in% "Primary solid Tumor")
            samplenr[snrcount, 2] <- sum(class %in% "Solid Tissue Normal")
        }
    }

    stop("tcga done")

}

## commands for hpc:

system("scp ~/Documents/mnem/src/mm.cpp euler:")
system("scp ~/Documents/mnem/R/mnems.r euler:")
system("scp ~/Documents/mnem/R/mnems_low.r euler:")
system("scp ~/Documents/nempi/R/nempi_main.r euler:")
system("scp ~/Documents/nempi/R/nempi_low.r euler:")

system("scp testing/nempi/other/nempi_loo.r euler:")

do="loo"

rm error.txt
rm output.txt

ram=8000 # nempi 8000, rf 64000, svm 20000, mf 8000, knn 8000, nn 20000
queue=4 # nempi 2min, rf 90min, svm 10min, mf 2min, knn 1min, nn 40min
bsruns=10
startrun=1
R="R/bin/R"

type="nn" # nempi, svm, rf, nn, knn, mf

if [ ${type} == 'rf' ] || [ ${type} == 'svm' ] || [ ${type} == 'nn' ]
then
maxpp=500000
else
maxpp=50000
fi
if [ ${type} == 'rf' ]
then
ram=64000
bsruns=1
fi
if [ ${type} == 'svm' ]
then
ram=20000
fi
if [ ${type} == 'nn' ]
then
ram=20000
bsruns=1
fi

bsub -M ${ram} -q normal.${queue}h -n 1 -e error.txt -o output.txt -R "rusage[mem=${ram}]" "${R} --vanilla --max-ppsize ${maxpp} --silent --no-save --args '${do}' '${type}' '${bsruns}' '${startrun}' < nempi_loo.r"

end=$(expr 91 + $bsruns )
end=$(expr $end / $bsruns )
## end=85

for (( i=2; i<=$end; i++ )); do
a=$(expr ${i} - 1 )
b=$(expr ${a} \* ${bsruns} )
startrun=$(expr ${b} + 1 )
echo $startrun
bsub -M ${ram} -q normal.${queue}h -n 1 -e error.txt -o output.txt -R "rusage[mem=${ram}]" "${R} --vanilla --max-ppsize ${maxpp} --silent --no-save --args '${do}' '${type}' '${bsruns}' '${startrun}' < nempi_loo.r"
done

for i in {2..10}; do
bsub -M ${ram} -q normal.${queue}h -n 1 -e error.txt -o output.txt -R "rusage[mem=${ram}]" "${R} --vanilla --max-ppsize ${maxpp} --silent --no-save --args '${do}' '${type}' '${bsruns}' '${startrun}' < nempi_loo.r"
done

## read in loo results:

methods <- c("nempi", "svm", "nn", "rf", "mf", "knn", "rand")
methlist <- list()
comparePre <- compareRec <- comparePrauc <- NULL
for (method in methods) {
    print(method)
    methlist[[method]] <- list()
    tp2 <- tn2 <- fp2 <- fn2 <- prauc2 <- numeric(91)
    for (file in list.files("~/Mount/Eulershare/nempi_bs")) {
        if (length(grep(paste0("^", method, "_"), file)) > 0 &
            length(grep(paste0("_loo\\."), file)) > 0) {
            load(paste0("~/Mount/Eulershare/nempi_bs/", file))
            tp2 <- tp2 + tp
            tn2 <- tn2 + tn
            fp2 <- fp2 + fp
            fn2 <- fn2 + fn
            prauc2 <- prauc2 + prauc
        }
    }
    methlist[[method]][["tp"]] <- tp2
    methlist[[method]][["tn"]] <- tn2
    methlist[[method]][["fp"]] <- fp2
    methlist[[method]][["fn"]] <- fn2
    methlist[[method]][["pre"]] <- tp2/(tp2+fp2)
    methlist[[method]][["rec"]] <- tp2/(tp2+fn2)
    methlist[[method]][["rec"]][which(is.na(methlist[[method]][["rec"]]))] <- 0
    methlist[[method]][["prauc"]] <- prauc2
    comparePre <- cbind(comparePre, methlist[[method]][["pre"]])
    compareRec <- cbind(compareRec, methlist[[method]][["rec"]])
    comparePrauc <- cbind(comparePrauc, methlist[[method]][["prauc"]])
}

setEPS()
postscript("temp.eps", height = 7, width = 10)
layout.mat <- matrix(c(rep(c(1,4,7),each=2),10,rep(c(2,5,8),each=2),11,rep(c(3,6,9),each=2),12),7)
par(mfrow=c(1,2))
col <- c("red", "blue", "darkgreen", "brown", "orange", "turquoise","grey")
mnem:::myboxplot(comparePre, dens = 0, col = col, xaxt = "n", main = "precision of mutation calls", ylab = "precision",box=1,border = col,medcol="black")
axis(1, 1:7, round(apply(comparePre,2,mean), 2))
mnem:::myboxplot(compareRec, dens = 0, col = col, xaxt = "n", main = "recall of mutation calls", ylab = "recall",box=1,border = col,medcol="black")
axis(1, 1:7, round(apply(compareRec,2,mean), 2))
dev.off()


setEPS()
postscript("temp.eps", height = 5, width = 5)
layout.mat <- matrix(c(rep(rep(1,each=3),3),2:4),4,byrow=1)
layout(layout.mat)
wtest <- numeric(ncol(comparePrauc))
for (i in 1:length(wtest)) {
    wtest[i] <- wilcox.test(comparePrauc[,1],comparePrauc[,i],alternative="greater")$p.value
}
col <- c("red", "blue", "darkgreen", "brown", "orange", "turquoise","grey")
mnem:::myboxplot(comparePrauc, dens = 0, col = col, xaxt = "n", main = "area under the precision-recall curve of mutation calls", ylab = "area under the precision-recall curve",box=1,border = col,medcol="black")
axis(1, 1:7, round(apply(comparePrauc,2,mean), 2))
par(mar=rep(0, 4))
plot.new()
legend("topleft",legend=c(expression(NEM~pi), "svm", "neural net"),col=c("red", "blue", "darkgreen"),fill=c("red", "blue", "darkgreen"),box.col = "transparent",cex = cex.leg,ncol = 1,border="transparent")
plot.new()
legend("topleft",legend=c("random forest", "missForest", "knn"),col=c("brown", "orange", "turquoise"),fill=c("brown", "orange", "turquoise"),box.col = "transparent",cex = cex.leg,ncol = 1,border="transparent")
plot.new()
legend("topleft",legend=c("random"),col=c("grey"),fill=c("grey"),box.col = "transparent",cex = cex.leg,ncol = 1,border="transparent")
dev.off()