# Random Forest classifier using RNA gene expression, mutation and both combined data

library(glmnet)
library(randomForest)
library(caret)
library(e1071)

source("plot_roc_curve.R")
options(warn=-1)


# Classification using RNA expression data
Classifier.RF.RNA <- function(input.RNA.file){

	rna.exp <- read.table(input.RNA.file, header = TRUE, row.names = 1)
	y <- rna.exp[,dim(rna.exp)[2]]
	x.in <- rna.exp[,c(1:(dim(rna.exp)[2]-1))]

	y.binary <- rep(0, length(y))
	y.binary[which(y=="YES")] <- 1

	prediction.values <- vector()

	folds <- cut(seq(1,nrow(x.in)),breaks=10,labels=FALSE)

	for(i in 1:10){
    		test.Indexes <- which(folds==i, arr.ind=TRUE)
    		print(test.Indexes)

    		trainData <- x.in[-test.Indexes, ]
    		train.y <-  y[-test.Indexes]
    		train.y.binary <- y.binary[-test.Indexes]
    		testData <- x.in[test.Indexes, ] 
    		test.y <- y[test.Indexes]
    		test.y.binary <- y.binary[test.Indexes]

    		cv.glmmod <- cv.glmnet(as.matrix(trainData), train.y.binary)
    		coef.result <- coef(cv.glmmod, s="lambda.min")
    		plot(cv.glmmod)

    		coef_inx <- summary(coef.result)[,1]-1
    		print(colnames(x.in)[coef_inx])

    		train.new <- data.frame(trainData[,coef_inx], train.y)
    		colnames(train.new)[dim(train.new)[2]] <- "CloneNumber"
    		test.new <- data.frame(testData[,coef_inx], test.y)
    		colnames(test.new)[dim(test.new)[2]] <- "CloneNumber"

    		# parameter tunning & training
    		nodesize.vals <- c(2,3,4,5)
    		ntree.vals <- c(200, 500, 1000, 2000)
    		tuning.results <- tune.randomForest(CloneNumber~.,data=train.new, mtry=3, 
			ntree=ntree.vals, nodesize=nodesize.vals)
    		clone_rf <- tuning.results$best.model
    		print(tuning.results)

    		importance <- varImp(clone_rf, scale=FALSE)
    		print(importance)

    		#print(clone_rf)

    		# testset prediction
    		clonePred.prob<-predict(clone_rf,newdata=test.new, type="prob")
    		values <- clonePred.prob[,2]
    		prediction.values <- c(prediction.values, values)

	}

	predictions <- prediction(prediction.values, y.binary)
	plot.roc.curve(predictions, title.text = "ROC Curve")

	predictions.round <- round(prediction.values)
	my_list <- list("pred" = as.factor(predictions.round), "target" = as.factor(y.binary) )
	return(my_list)
}
  
# Classification using mutation data
Classifier.RF.Mutation <- function(input.Mutation.file){

	dat.mut <- read.table(input.Mutation.file, header = TRUE, row.names = 1)
	y <- dat.mut[,dim(dat.mut)[2]]
	x.in <- dat.mut[,c(1:(dim(dat.mut)[2]-1))]

	y.binary <- rep(0, length(y))
	y.binary[which(y=="YES")] <- 1

	prediction.values <- vector()

	folds <- cut(seq(1,nrow(x.in)),breaks=10,labels=FALSE)

	for(i in 1:10){
    		test.Indexes <- which(folds==i, arr.ind=TRUE)
    		print(test.Indexes)

    		trainData <- x.in[-test.Indexes, ]
    		train.y <-  y[-test.Indexes]
    		train.y.binary <- y.binary[-test.Indexes]
    		testData <- x.in[test.Indexes, ] 
    		test.y <- y[test.Indexes]
    		test.y.binary <- y.binary[test.Indexes]

    		cv.glmmod <- cv.glmnet(as.matrix(trainData), train.y.binary)
    		coef.result <- coef(cv.glmmod, s="lambda.min")
    		plot(cv.glmmod)

    		coef_inx <- summary(coef.result)[,1]-1
    		print(colnames(x.in)[coef_inx])

    		train.new <- data.frame(trainData[,coef_inx], train.y)
    		colnames(train.new)[dim(train.new)[2]] <- "CloneNumber"
    		test.new <- data.frame(testData[,coef_inx], test.y)
    		colnames(test.new)[dim(test.new)[2]] <- "CloneNumber"

    		# parameter tunning & training
    		nodesize.vals <- c(2,3,4,5)
    		ntree.vals <- c(200, 500, 1000, 2000)
    		tuning.results <- tune.randomForest(CloneNumber~.,data=train.new, mtry=3, 
                         ntree=ntree.vals, nodesize=nodesize.vals)
    		clone_rf <- tuning.results$best.model
    		print(tuning.results)

    		importance <- varImp(clone_rf, scale=FALSE)
    		print(importance)

    		#print(clone_rf)

    		# testset prediction
    		clonePred.prob<-predict(clone_rf,newdata=test.new, type="prob")
    		values <- clonePred.prob[,2]
    		prediction.values <- c(prediction.values, values)

	}

	predictions <- prediction(prediction.values, y.binary)
	plot.roc.curve(predictions, title.text = "ROC Curve")

	predictions.round <- round(prediction.values)
	my_list <- list("pred" = as.factor(predictions.round), "target" = as.factor(y.binary) )
	return(my_list)
}

# Classification using combined data
Classifier.RF.Combined <- function(input.RNA.file, input.Mutation.file){

	rna.exp <- read.table(input.RNA.file, header = TRUE, row.names = 1)
	x.in.rna <- rna.exp[,c(1:(dim(rna.exp)[2]-1))]
	y <- rna.exp[,dim(rna.exp)[2]]

	dat.mut <- read.table(input.Mutation.file, header = TRUE, row.names = 1)
	x.in.mut <- dat.mut[,c(1:(dim(dat.mut)[2]-1))]

	y.binary <- rep(0, length(y))
	y.binary[which(y=="YES")] <- 1

	prediction.values <- vector()

	folds <- cut(seq(1,nrow(x.in.rna)),breaks=10,labels=FALSE)

	for(i in 1:10){
    		test.Indexes <- which(folds==i, arr.ind=TRUE)
    		print(test.Indexes)

    		trainData.rna <- x.in.rna[-test.Indexes, ]
    		testData.rna <- x.in.rna[test.Indexes, ] 
    		trainData.mut <- x.in.mut[-test.Indexes, ]
    		testData.mut <- x.in.mut[test.Indexes, ] 

    		train.y <-  y[-test.Indexes]
    		train.y.binary <- y.binary[-test.Indexes]
    		test.y <- y[test.Indexes]
    		test.y.binary <- y.binary[test.Indexes]

    		# Lasso for RNA
    		cv.glmmod.rna <- cv.glmnet(as.matrix(trainData.rna), train.y.binary)
    		coef.result.rna <- coef(cv.glmmod.rna, s="lambda.min")
    		plot(cv.glmmod.rna)

    		coef_inx.rna <- summary(coef.result.rna)[,1]-1
    		print(colnames(x.in.rna)[coef_inx.rna])

    		# Lasso for mutation
    		cv.glmmod.mut <- cv.glmnet(as.matrix(trainData.mut), train.y.binary)
    		coef.result.mut <- coef(cv.glmmod.mut, s="lambda.min")
    		plot(cv.glmmod.mut)

    		coef_inx.mut <- summary(coef.result.mut)[,1]-1
    		print(colnames(x.in.mut)[coef_inx.mut])

    		# Combined dataset
    		train.new <- data.frame(trainData.rna[,coef_inx.rna], trainData.mut[,coef_inx.mut], train.y)
    		colnames(train.new)[dim(train.new)[2]] <- "CloneNumber"
    		test.new <- data.frame(testData.rna[,coef_inx.rna], testData.mut[,coef_inx.mut], test.y)
    		colnames(test.new)[dim(test.new)[2]] <- "CloneNumber"

    		# parameter tunning & training
    		nodesize.vals <- c(2,3,4,5)
    		ntree.vals <- c(200, 500, 1000, 2000)
    		tuning.results <- tune.randomForest(CloneNumber~.,data=train.new, mtry=3, ntree=ntree.vals, 
                           nodesize=nodesize.vals)
    		clone_rf <- tuning.results$best.model
    		print(tuning.results)

    		importance <- varImp(clone_rf, scale=FALSE)
    		print(importance)

    		#print(clone_rf)

    		clonePred.prob<-predict(clone_rf,newdata=test.new, type="prob")
    		values <- clonePred.prob[,2]
    		prediction.values <- c(prediction.values, values)

	}

	predictions <- prediction(prediction.values, y.binary)
	plot.roc.curve(predictions, title.text = "ROC Curve")

	predictions.round <- round(prediction.values)
	my_list <- list("pred" = as.factor(predictions.round), "target" = as.factor(y.binary) )
	return(my_list)
}