Step 3：Machine Learning Screening for Hub Genes
setwd("-")
library(tidyverse)
library(glmnet)
source('msvmRFE.R')
library(e1071)
library(caret)
library(bRacatus)
library(randomForest)
gene<-read.table('gene.txt',sep='\t',header=T)
names(gene)<-'gene'
load("SD_exp_combined-human.Rdata")
exp2=exp2[,-42]
exp <- merge(gene,exp2,by.x = 1,by.y = 0)
rownames(exp)=exp[,1]
exp=exp[,-1]
exp=as.data.frame(t(exp))

set.seed(12345) 
x<-as.matrix(exp)
head(group_list)
y<-ifelse(group_list=='Control',0,1)
fit<-glmnet(x,y,family = 'binomial',alpha = 1,lambda = NULL)
pdf('lasso.pdf',width = 8,height = 7)
plot(fit,xvar = 'lambda',label=T)
dev.off()
cvfit<-cv.glmnet(x,y,family='binomial',alpha = 1,type.measure = 'deviance',nfolds = 5)
pdf('cvfit.pdf',width = 8,height = 7)
plot(cvfit)
dev.off()
best_lambda<-cvfit$lambda.min
lasso_coefs<-coef(cvfit,s='lambda.min')
lasso_features<-lasso_coefs@Dimnames[[1]][which(lasso_coefs!=0)]
lasso_features=as.data.frame(lasso_features)
lasso_features<-lasso_features[-1,]
write.csv(lasso_features,'feature_lasso.csv')
train<-cbind(y,exp)

colnames(train)[1]<-'group'
input<-train
nfold<-5   
nrows<-nrow(input)
folds<-rep(1:nfold,len=nrows)[sample(nrows)]
folds<-lapply(1:nfold,function(x) which(folds==x))
results<-lapply(folds,svmRFE.wrap,input,k=5,halve.above=20)
top.features<-WriteFeatures(results,input,save=F)
write.csv(top.features,'feature_svm.csv')

FeatSweep<-lapply(1:5,FeatSweep.wrap,results,input)
no.info<-min(prop.table(table(input[,1])))
errors<-sapply(FeatSweep,function(x) ifelse(is.null(x),NA,x$error))
pdf('svm-error.pdf',width = 6,height = 6)
PlotErrors(errors,no.info = no.info)
dev.off()
Plotaccuracy<-function(accuracies,no.info=0.5,ylim=range(accuracies),
                       xlab='Number of Features',ylab='5x CV Accuracy'){
  AddLine<-function(x,col='black'){
    max_index<-which.max(x)
    lines(x=max_index,y=x[max_index],col=col,type='p')
    points(x=max_index,y=x[max_index],col='red')
    text(x=max_index,y=x[max_index],labels=paste(max_index,'-',format(x[max_index],digits=3)),pos=1,
         col='red',cex=0.75)
  }
  plot(x=1:length(accuracies),y=accuracies,type='l',ylim=ylim,xlab=xlab,ylab=ylab)
  AddLine(accuracies)
  abline(h=no.info,lty=3)
}

pdf('svm-accuracy.pdf',width = 8,height = 8)
Plotaccuracy(1-errors,no.info = no.info)
dev.off()
min_error_index<-which.min(errors)
write.csv(SVMRFEgenes,file='SVMRFEgenes.csv')


library(ggvenn)
library(ggplot2)
library(dplyr)
Lasso<-read.csv('feature_lasso.csv',sep=',',header = T)
SVMRFE<-read.csv('SVMRFEgenes.csv',sep=',',header = T)
datalist<-list('Lasso'=Lasso$x,
               'SVM-RFE'=SVMRFE$x)
opar<-par(family='Roboto Condensed')
biocolor<-c('blue','yellow')
p=ggvenn(datalist,
         fill_color = biocolor,
         fill_alpha = 0.5,
         stroke_linetype = 'longdash',
         set_name_size = 3,
         text_size = 4)
p
pdf(file='ML.pdf',width = 8,height = 6)
print(p)
dev.off()


gene_list<-list(Lasso$x,SVMRFE$x)
common_gene<-Reduce(intersect,gene_list)
write.csv(common_gene,'ML.csv')