R/GenoNet.R

In Ionita-Laza-lab/GenoNet: A Semi-Supervised Approach for Predicting Cell-Type Specific Functional Consequences of Non-Coding Variation

GenoNet<-function(Y,X,beta.unlabeled,nfolds=10){

  Y<-as.matrix(Y);X<-as.matrix(X);beta.unlabeled<-as.matrix(beta.unlabeled)
  X<-cbind(1,X) #add intercept term
  #shuffle the samples for cross-validation
  index.shuffle<-sample(1:length(Y))
  Y<-as.matrix(Y[index.shuffle,]);X<-as.matrix(X[index.shuffle,])

  lambda.opt<-sd(Y)*sqrt(2*log(ncol(X)-1)/nrow(X))
  list.ENET<-list()
  for(i in 1:nfolds){
    valid.index<-(i-1)*floor(length(Y)/nfolds)+1:floor(length(Y)/nfolds)
    valid.X<-X[valid.index,];valid.Y<-as.matrix(Y[valid.index])
    beta0.ENET<-as.matrix(coef(glmnet(X[-valid.index,-1],Y[-valid.index],family='binomial',standardize=F,alpha=0.5,lambda=lambda.opt)))
    list.ENET[[i]]<-beta0.ENET
  }

  if(length(Y)<=1000){
    Cross.COR<-function(phi){
      cross.Y<-c();cross.score<-c()
      for(i in 1:nfolds){
        valid.index<-(i-1)*floor(length(Y)/nfolds)+1:floor(length(Y)/nfolds)
        valid.X<-X[valid.index,];valid.Y<-as.matrix(Y[valid.index])
        beta0.ENET<-list.ENET[[i]]
        mu<-valid.X%*%beta0.ENET
        ENET.valid<-exp(mu)/(1+exp(mu))
        L1kg.valid<-valid.X%*%beta.unlabeled
        L1kg.valid<-L1kg.valid*(L1kg.valid<=1 & L1kg.valid>=0)+(L1kg.valid>1)
        valid.score<-(1-phi)*ENET.valid+phi*L1kg.valid
        cross.Y<-c(cross.Y,valid.Y);cross.score<-c(cross.score,valid.score)
      }
      cross.COR<--abs(cor(cross.Y,cross.score))
      return(cross.COR)
    }

    fit.phi<-optim(par=1,Cross.COR,method='L-BFGS-B',lower=0,upper=1)
    phi<-fit.phi$par
  }else{
    Cross.AUC<-function(phi){
      cross.Y<-c();cross.score<-c()
      for(i in 1:nfolds){
        valid.index<-(i-1)*floor(length(Y)/nfolds)+1:floor(length(Y)/nfolds)
        valid.X<-X[valid.index,];valid.Y<-as.matrix(Y[valid.index])
        beta0.ENET<-list.ENET[[i]]
        mu<-valid.X%*%beta0.ENET
        ENET.valid<-exp(mu)/(1+exp(mu))
        L1kg.valid<-valid.X%*%beta.unlabeled
        L1kg.valid<-L1kg.valid*(L1kg.valid<=1 & L1kg.valid>=0)+(L1kg.valid>1)
        valid.score<-(1-phi)*ENET.valid+phi*L1kg.valid
        cross.Y<-c(cross.Y,valid.Y);cross.score<-c(cross.score,valid.score)
      }
      cross.AUC<--abs(Get.AUC(cross.score,cross.Y))
      return(cross.AUC)
    }

    fit.phi<-c(Cross.AUC(0),Cross.AUC(1))
    phi<-c(0,1)[which.min(fit.phi)]
  }

  lambda.opt<-sd(Y)*sqrt(2*log(ncol(X)-1)/nrow(X))
  beta.labeled<-as.matrix(coef(glmnet(X[,-1],Y,family='binomial',standardize=F,alpha=0.5,lambda=lambda.opt)))
  return(list(LogitBeta.labeled=beta.labeled,LinearBeta.unlabeled=beta.unlabeled,phi=phi))
}

GenoNet.predict<-function(X,GenoNet.fit){
  X<-as.matrix(X)
  X<-cbind(1,X)

  mu<-X%*%GenoNet.fit$LogitBeta.labeled
  labeled.predict<-exp(mu)/(1+exp(mu))

  unlabeled.predict<-X%*%GenoNet.fit$LogitBeta.labeled
  unlabeled.predict<-unlabeled.predict*(unlabeled.predict<=1 & unlabeled.predict>=0)+(unlabeled.predict>1)

  Y.predict<-(1-GenoNet.fit$phi)*labeled.predict+GenoNet.fit$phi*unlabeled.predict

  return(Y.predict=Y.predict)
}

Get.AUC<-function(x,y){
  df <- data.frame(x=x, y=y)
  utest <- wilcox.test(x ~ y, data=df)#utest$statistic / prod(table(df$y))
  AUC<-max(utest$statistic / prod(table(df$y)),1 - utest$statistic / prod(table(df$y)))
  return(AUC)
}