E2Predictor:

library(E2Predictor)
library(gtools)
library(ggplot2)

### Configuration #########################
initConfiguration(netMHCpath = "/Users/napedro/external_sources/netMHC-4.0/netMHC",
                  working.path = "/Users/napedro/CloudStation")

dbfile <- file.path(E2Predictor.Config$working.path, "databases", "canonical", "uniprot-human.tab.gz")

class_category = "has_ligands"

cell.lines <- names(E2Predictor.Config$hla_alleles)
conditions <- c("Control", "Bortezomib", "DMOG", "IFNg", "Rapamycin")
cell.line.training <- cell.lines[1]
cell.line.testing <- cell.lines[2]

is_perturbation <- F
condition <- conditions[5]
addgrepcond <- ""  # " II "

analysis_folder <- file.path(E2Predictor.Config$working.path, "data_analyses", "machine_learning")

requireAllele = TRUE
factor.training <- 0.6
hidden.layers <- c(6, 2)

set.seed(1234567890)
###########################################

## load data: this loads tables from the folder "tables.for.modeling". They include protein values, RPKM, alleles, and GO enriched terms
omics.cellline1 <- load_data_modeling(cell.line = cell.line.training, is_perturbation = F, requireAllele = requireAllele)
omics.cellline2 <- load_data_modeling(cell.line = cell.line.testing, is_perturbation = F, requireAllele = requireAllele)
omics.combined <- rbind(omics.cellline1, omics.cellline2)
rm(omics.cellline1, omics.cellline2)

## additional knowledge databases

# crawler includes several predictions for cell localization and presence of helices
crawler <- read_file_crawler( file.path(E2Predictor.Config$working.path, "predictions.crowler/Uniprot_Predictions_Localization_Crawler.xlsx") )
crawler.preML <- prelim_ML(df = crawler, dim.reduction.factor = 0.6, save.to = file.path(analysis_folder, "crawler"))
crawler.red <- crawler[, !names(crawler) %in% crawler.preML$highly.correlated.vars]

# SwissProt database, including many info about each protein. Here I filtered it to only 5
# features (Length, Mass, Transmembrane, Glycosylation and Lipidation), but you may include more
swissprot <- readr::read_tsv(file = dbfile)
swissprot <- swissprot %>%
    filter(Status == "reviewed") %>%
    rename(protein_entry = Entry) %>%
    select(protein_entry, Length, Mass, Transmembrane, Glycosylation, Lipidation) %>%
    mutate( Transmembrane = ifelse(is.na(Transmembrane), FALSE, TRUE)
            ,Glycosylation = ifelse(is.na(Glycosylation), FALSE, TRUE)
            ,Lipidation = ifelse(is.na(Lipidation), FALSE, TRUE))


# Table generated by Sebastian with netctl predictions
netctl <- read.table(file.path(E2Predictor.Config$working.path, "predictions", "netctl.tsv"), sep = "\t", header = T)
# netctl.preML <- prelim_ML(netctl, dim.reduction.factor = 0.6, save.to = file.path(analysis_folder, "netctl_table"))

# Merge several tables. Remember that here you can also add more cell lines, etc...
omics.combined <- merge(omics.combined, crawler.red, by="protein_entry")
omics.combined <- merge(omics.combined, swissprot, by="protein_entry")
omics.combined <- merge(omics.combined, netctl, by.x="Entry name", by.y = "Entry.name")

# prepare data modeling

## Defined variables for preparing data modeling. I guess it would be enough to have just one definition for these variables,
## so that all numeric (and boolean) variables are normalized. After wards you may train over a subset of variables by using
## a different formula
df1.vars <- c("has_ligands", "PROTEIN_QUANTITY", "RPKM.avg",  grep("\\[GO.*", names(omics.combined), value = T)) #
df2.vars <- c(df1.vars, "Length", "Mass", "Transmembrane", "Glycosylation", "Lipidation", names(crawler.red)[2:length(crawler.red)])
df.vars <- c(df2.vars, names(netctl)[2:length(names(netctl))])
df.old.vars <- c(df1.vars, names(netctl)[2:length(names(netctl))])

## Prepare the data for modeling, by using only the variables defined above (I still keep these three variable subsets, which is useless...)
omics.combined.datamodel <- prepare_data_modeling(omics.combined, useVariables = df.vars)
omics.combined.datamodel.old <- prepare_data_modeling(omics.combined, useVariables = df.old.vars)

## test / training datasets. This returns an object containing two data frames: trainset and testset
omics.combined_tr <- get_training_sample(omics.combined.datamodel, class_category = class_category, factor.training = 1, grep_rowname = "^JY")
omics.combined_tr.old <- get_training_sample(omics.combined.datamodel.old, class_category = class_category, factor.training = 1, grep_rowname = "^JY")
omics.combined_tr_inverse <- get_training_sample(omics.combined.datamodel, class_category = class_category, factor.training = 1, grep_rowname = "^LCLC")

## Preliminary Machine Learning (var. correlation & PCA)
omics.combined.preML <- prelim_ML(df= omics.combined.datamodel, dim.reduction.factor = 0.7, save.to = file.path(analysis_folder, "test"))

## Machine Learning

f.simple <- as.formula(paste("has_ligands ~ PROTEIN_QUANTITY + RPKM.avg"))
f.all <- as.formula(paste(class_category, "~", paste(names(omics.combined_tr$trainset)[!names(omics.combined_tr$trainset) %in% class_category], collapse = " + ")))
f.all.old <- as.formula(paste(class_category, "~", paste(names(omics.combined_tr.old$trainset)[!names(omics.combined_tr.old$trainset) %in% class_category], collapse = " + ")))

decTreeFit.all <- rpart::rpart(f.all, data = omics.combined_tr$trainset, method = "anova")
decTreeFit.all.old <- rpart::rpart(f.all.old, data = omics.combined_tr$trainset, method = "anova")

rpart::printcp(decTreeFit.all)
rpart::plotcp(decTreeFit.all)
plot(decTreeFit.all, uniform=TRUE)
summary(decTreeFit.all)
summary(decTreeFit.all.old)

decTreeFit.all$variable.importance
decTreeFit.all.old$variable.importance

## Evaluate Machine Learning. It saves some results (ROC curves, etc) at a user defined folder
path.decTreeFit.all = file.path(analysis_folder, "decTree.all")
path.ANN.Comb = file.path(analysis_folder, "combined.JY.LCLC_ANN")
path.ANN.Optimization.base = file.path(analysis_folder, "ANN.JYandLCLC.optimization")

features <- c("RPKM.avg"
              ,"PROTEIN_QUANTITY"
              ,"Glycosylation"
              ,"TARGETP.Location.S"
              ,"SUBLOC.Location.Extracellular"
              ,"GO0004252"
              ,"GO0005576"
              ,"GO0004930"
              ,"tap_1.0"
              ,"tap_0.5"
              ,"cleavage_0.5"
              ,"cleavage_0.7"
              ,"cleavage_0.9"
#              ,"TMPRED.Helices"
              ,"binder_rel_A2"
            )

feat.combinations <- combinations(n = length(features), r=length(features) -1 , v=features, set=TRUE, repeats.allowed = FALSE)

run.ANN <- function(feats, traintest, folder.name){

    f <- as.formula( paste(class_category, "~" ,  paste(feats, collapse = "+")) )
    print(paste0("Running ANN: ", n.comb))
    print(f)
    save.to <- file.path(path.ANN.Optimization.base, folder.name, paste("comb", n.comb, sep = "_"))
    mkdir(save.to)
    sink(file = file.path(save.to, "features.txt"))
        print(f)
    sink()

    annFit <- neuralnet::neuralnet(f, data=traintest$trainset, hidden=hidden.layers, linear.output=F, stepmax = 1e+06)
    annFit.ROC <- evaluate_ML(predictor=annFit, testset=traintest$testset, class_category = class_category, save.to = save.to)

    n.comb <<- n.comb + 1

    return(list( predictor = annFit, roc = annFit.ROC) )
}

run.ANN.nnet <- function(feats, traintest, folder.name){

    f <- as.formula( paste(class_category, "~" ,  paste(feats, collapse = "+")) )
    print(paste0("Running ANN: ", n.comb))
    print(f)
    save.to <- file.path(path.ANN.Optimization.base, folder.name, paste("comb", n.comb, sep = "_"))
    mkdir(save.to)
    sink(file = file.path(save.to, "features.txt"))
    print(f)
    sink()

    annFit <- nnet::nnet(f, data=traintest$trainset, size = 10, linout=F, skip=T, maxit=10000)
    if(is.null(annFit)) return(NA)
    annFit.ROC <- evaluate_ML(predictor=annFit, testset=traintest$testset, class_category = class_category, save.to = save.to)


    p1<-NeuralNetTools::lekprofile(annFit)

    p1_tr <-p1 +
        theme_bw() +
        scale_colour_brewer(palette="PuBu") +
        scale_linetype_manual(values=rep('dashed',6)) +
        scale_size_manual(values=rep(1,6))

    ggsave(file.path(save.to, "lek.plot.pdf"), plot = p1_tr, width = 30, height = 7)

    n.comb <<- n.comb + 1

    return(list( predictor = annFit, roc = annFit.ROC) )
}


# n.comb <- 1
# permutations.training.JY <- apply(feat.combinations, 1, run.ANN, omics.combined_tr, "permutations.training.JY")
#
# taps.out.features <- features[grep("^tap", features, invert = TRUE)]
# n.comb <- 1
# permutations.training.JY <- run.ANN(feats = taps.out.features, omics.combined_tr, "no.taps.JY")
#
# cleavages.out.features <- features[grep("^cleav", features, invert = TRUE)]
# n.comb <- 1
# permutations.training.JY <- run.ANN(feats = cleavages.out.features, omics.combined_tr, "no.cleavages.JY")
#
# GO.out.features <- features[grep("^GO", features, invert = TRUE)]
# n.comb <- 1
# permutations.training.JY <- run.ANN(feats = GO.out.features, omics.combined_tr, "no.GOs.JY")
#
# localiz.out.features <- features[grep("Location", features, invert = TRUE)]
# n.comb <- 1
# permutations.training.JY <- run.ANN(feats = localiz.out.features, omics.combined_tr, "no.TARGETPandSUBLOC.JY")



# n.comb <- 1
# permutations.training.LCLC <- apply(feat.combinations, 1, run.ANN, omics.combined_tr_inverse, "permutations.training.LCLC")
#
# taps.out.features <- features[grep("^tap", features, invert = TRUE)]
# n.comb <- 1
# permutations.training.LCLC <- run.ANN(feats = taps.out.features, omics.combined_tr_inverse, "no.taps.LCLC")
#
# cleavages.out.features <- features[grep("^cleav", features, invert = TRUE)]
# n.comb <- 1
# permutations.training.LCLC <- run.ANN(feats = cleavages.out.features, omics.combined_tr_inverse, "no.cleavages.LCLC")
#
# GO.out.features <- features[grep("^GO", features, invert = TRUE)]
# n.comb <- 1
# permutations.training.LCLC <- run.ANN(feats = GO.out.features, omics.combined_tr_inverse, "no.GOs.LCLC")
#
# localiz.out.features <- features[grep("Location", features, invert = TRUE)]
# n.comb <- 1
# permutations.training.LCLC <- run.ANN(feats = localiz.out.features, omics.combined_tr_inverse, "no.TARGETPandSUBLOC.LCLC")



n.comb <- 1
permutations.training.nnet.JY <- apply(feat.combinations, 1, run.ANN.nnet, omics.combined_tr, "permutations.training.nnet.skipTrue.JY")


library(clusterGeneration)

seed.val<-2
set.seed(seed.val)

num.vars<-8
num.obs<-1000

#input variables
cov.mat<-genPositiveDefMat(num.vars,covMethod=c("unifcorrmat"))$Sigma
rand.vars<-mvrnorm(num.obs,rep(0,num.vars),Sigma=cov.mat)

#output variables
parms<-runif(num.vars,-10,10)
y1<-rand.vars %*% matrix(parms) + rnorm(num.obs,sd=20)
parms2<-runif(num.vars,-10,10)
y2<-rand.vars %*% matrix(parms2) + rnorm(num.obs,sd=20)

#final datasets
rand.vars<-data.frame(rand.vars)
resp<-data.frame(y1,y2)
names(resp)<-c('Y1','Y2')
dat.in<-data.frame(resp,rand.vars)

#import the function from Github
library(devtools)
source_url('https://gist.githubusercontent.com/fawda123/7471137/raw/466c1474d0a505ff044412703516c34f1a4684a5/nnet_plot_update.r')
library(NeuralNetTools)
library(RSNNS)

#neural net with two hidden layers, 9 and 2 nodes in each
trainset <- omics.combined_tr$trainset[, features]
traintarget <- omics.combined_tr$trainset[, (names(omics.combined_tr$trainset) %in% class_category)]
testset <- omics.combined_tr$testset[, features]
testtarget <- omics.combined_tr$testset[, (names(omics.combined_tr$testset) %in% class_category)]

traintarget <- factor(traintarget,
       levels=c(TRUE,FALSE),
       labels=c(TRUE,FALSE))

mod<-RSNNS::mlp(trainset, traintarget, size=c(9,2),linOut=T, inputsTest = testset, targetsTest = testtarget)
mod.roc <- evaluate_ML(mod, testset, class_category = class_category, save.to = file.path(path.ANN.Optimization.base, "test_mlp"))


# pruned model using code from RSSNS pruning demo
pruneFuncParams <- list(
    max_pr_error_increase = 10.0,
    pr_accepted_error = 1.0,
    no_of_pr_retrain_cycles = 1000,
    min_error_to_stop = 0.01,
    init_matrix_value = 1e-6,
    input_pruning = TRUE,
    hidden_pruning = TRUE)
mod <- RSNNS::mlp(trainset, traintarget, size = c(10,4), learnFunc="BackpropMomentum", maxit= 10000)
        #pruneFunc = "OptimalBrainSurgeon", pruneFuncParams = pruneFuncParams)
mod.roc <- evaluate_ML(mod, testset, class_category = class_category,evaluation.features = features, save.to = file.path(path.ANN.Optimization.base, "test_mlp_BackpropMomentum"))


## caret ANN

all.feats <- names(omics.combined_tr$trainset)
mod.caret.nnet.20 <- caret::train(#trainset,
                                  omics.combined_tr$trainset[, !(names(omics.combined_tr$trainset) %in% class_category)],
                                     traintarget,
                                     #preProcess="pca",
                                     metric="Accuracy",
                                     maximize=T,
                                     method="nnet", maxit=10000, skip=TRUE # , size=20
                                )




mod.caret.nnet.roc <- evaluate_ML(predictor=mod.caret.nnet$finalModel,
                                  testset=testset,
                                  class_category=class_category,
                                  evaluation.features=features,
                                  save.to=file.path(path.ANN.Optimization.base, "test_mlp_caret_ann_allFeats_skipTrue"))





#dput(omics.combined_tr$trainset)
omics.combined_tr$trainset[,class_category] <- as.factor(omics.combined_tr$trainset[,class_category])
mod.caret.neuralnet.20.5.2 <- caret::train(
    form=f.all,
    data=omics.combined_tr$trainset,
    preProcess="pca",
    metric="Accuracy",
    maximize=T,
    method="neuralnet", layer1=20, layer2=5, layer3=2
)