##ssGSEA
library(GSVA)
library(GSVAdata)
geneSets <- getGmt("D:/GSEA/h.all.v6.2.symbols.gmt")
topMatrixGSVA <- gsva(as.matrix(rsem), geneSets, min.sz=10, max.sz=999999, abs.ranking=FALSE, verbose=TRUE)
topMatrixGSVA <- data.frame(topMatrixGSVA)
names(topMatrixGSVA) <- str_sub(names(topMatrixGSVA), 1, 15)
topMatrixGSVA <- data.frame(t(topMatrixGSVA))
rownames(topMatrixGSVA) <- str_replace_all(rownames(topMatrixGSVA), '[.]', '-')
topMatrixGSVA[1:3, 1:3]
names(topMatrixGSVA) <- str_split(names(topMatrixGSVA), 'HALLMARK_') %>%
  map(., function(x) x[[2]]) %>%
  unlist 
tem_topMatrixGSVA <- topMatrixGSVA
colnames <- colnames(topMatrixGSVA)
all(row.names(target) == rownames(topMatrixGSVA))
data <- data.frame(colnames)
for (j in 1:length(colnames)){
  test <- cor.test(as.numeric(topMatrixGSVA[,j]),target$MS4A15,type="spearman")
  data[j,2] <- test$estimate                                            
  data[j,3] <- test$p.value
  names(data) <- c("symbol","correlation","pvalue")
}
data <- data[order(data$correlation, decreasing = TRUE), ]
data %>% 
  filter(pvalue <0.01) %>% 
  ggplot(aes(correlation,forcats::fct_reorder(symbol,correlation))) +
  geom_segment(aes(xend=0,yend=symbol)) +
  geom_point(aes(col=pvalue,size=abs(correlation))) +
  scale_colour_gradientn(colours=c("#7fc97f","#984ea3")) +
  scale_size_continuous(range =c(2,8))  +
  theme_minimal() +
  ylab(NULL) + small
ggsave('ssGSEA.pdf', height = 9, width = 9)


##WGCNA
library(WGCNA)
Args <- commandArgs()
fFolder=as.character(Args[6])			
ePath=as.character(Args[7])			  	
sPath=as.character(Args[8])				  
cutSd=as.numeric(Args[9])			  	
minModuleSize=as.numeric(Args[10]) 
type=as.character(Args[11])        
BlockSize=as.numeric(Args[12])      
mergeCutHeight=as.numeric(Args[13])	
cutR=as.numeric(Args[14])			   
deepSplit=as.numeric(Args[15])     
cutNet=as.numeric(Args[16])         
cutHight=as.numeric(Args[17])        
samp=as.numeric(Args[18])  

#fFolder="D:\\test\\WGCNA_Test\\WGCNA6"
#ePath="success_top50var.symbol.exp.txt"
#sPath="success_SampleInfo_input.txt"
#cutSd=0
#minModuleSize=30
#type="unsigned"
#BlockSize=5000
#mergeCutHeight=0.25
#cutR=0.9
#deepSplit=2
#cutNet=0.02
#cutHight=-1
#samp=1


#cutHight=0;
#BlockSize=6000
#minModuleSize=30
#mergeCutHeight=0.25
#cutR=0.9
#cutNet=0.02
#cutSd=1.2
cex1=0.9
#ePath='206merge.FPKM.txt'
#sPath='206sample-group.txt'
#type='unsigned'
#deepSplit=4
#samp=1
files=c()
setwd(fFolder)
if(!require(WGCNA)){
  source("http://bioconductor.org/biocLite.R")
  options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
  biocLite('WGCNA')
  biocLite('rlang')
  biocLite('stringi')
  biocLite('reshape2')
  require(WGCNA)
}

corType='pearson'

print(ePath)
print(sPath)
#getwd()
expro=read.csv(ePath,sep = '\t',row.names = 1,stringsAsFactors = F,	check.names=F)
print(head(expro))
nclinical=0
if(sPath!=""&&sPath!=' '&&sPath!='0'){
  expro.sample=read.csv(sPath,sep = '\t',row.names = 1,stringsAsFactors = F, check.names=F)
  expro.sample=as.data.frame(expro.sample)
  nclinical=ncol(expro.sample)
}
print(sPath)
OutSampleCluster=function(sampleTree,expro.sample){
  if(nclinical>0){
    sColors=apply(expro.sample,2,function(x){
      cx=x
      ux=sort(unique(x))
      colors=colorRampPalette(c("green", "orange","red"))(length(ux))
      for(i in 1:length(ux)){
        cx[which(x==ux[i])]=colors[i]
      }
      cx[which(is.na(x))]="white"
      return(cx)
    })
    inds=c()
    for(i in 1:ncol(expro.sample)){
      if(class(expro.sample[,i])=='integer'||class(expro.sample[,i])=='numeric'){
        inds=c(inds,i)
      }
    }
    if(length(inds)==0){
      cols3=colorRampPalette(c("green", "orange","red"))(3)
      cols4=colorRampPalette(c("green", "orange","red"))(4)
      sColors=data.frame(Cluster2=ifelse(cutree(sampleTree,2)==1,'red','blue'),
                         Cluster3=cols3[cutree(sampleTree,3)],
                         Cluster4=cols4[cutree(sampleTree,4)]
      )
      expro.sample1=data.frame(Cluster2=cutree(sampleTree,2)
                               ,Cluster3=cutree(sampleTree,3)
                               ,Cluster4=cutree(sampleTree,4))
    }else{
      expro.sample1=as.data.frame(expro.sample[,inds])
      sColors=as.data.frame(sColors[,inds])
      colnames(sColors)=colnames(expro.sample)[inds]
      colnames(expro.sample1)=colnames(expro.sample)[inds]
    }
  }else{
    cols3=colorRampPalette(c("green", "orange","red"))(3)
    cols4=colorRampPalette(c("green", "orange","red"))(4)
    sColors=data.frame(Cluster2=ifelse(cutree(sampleTree,2)==1,'red','blue'),
                       Cluster3=cols3[cutree(sampleTree,3)],
                       Cluster4=cols4[cutree(sampleTree,4)]
    )
    expro.sample1=data.frame(Cluster2=cutree(sampleTree,2)
                             ,Cluster3=cutree(sampleTree,3)
                             ,Cluster4=cutree(sampleTree,4))
  }
  #sColors=NULL
  jpeg('SampleClusterAndColors.jpg',width = 800,height = 600)
  if(nrow(sColors)<50){
    plotDendroAndColors(sampleTree
                      ,sColors
                      ,colnames(sColors)
                      ,dendroLabels = NULL
                      ,addGuide = F
                      ,hang = 0.03
                      ,guideHang = 0.05)
  }else{
    plotDendroAndColors(sampleTree
                        ,sColors
                        ,colnames(sColors)
                        ,dendroLabels = F
                        ,addGuide = T
                        ,hang = 0.03
                        ,guideHang = 0.05)
  }
  dev.off()
  
  pdf('SampleClusterAndColors.pdf',width = 8,height = 6)
  if(nrow(sColors)<50){
    plotDendroAndColors(sampleTree
                        ,sColors
                        ,colnames(sColors)
                        ,dendroLabels = NULL
                        ,addGuide = F
                        ,hang = 0.03
                        ,guideHang = 0.05)
  }else{
    plotDendroAndColors(sampleTree
                        ,sColors
                        ,colnames(sColors)
                        ,dendroLabels = F
                        ,addGuide = T
                        ,hang = 0.03
                        ,guideHang = 0.05)
  }
  dev.off()
  return(expro.sample1)
}
checkClinicalData=function(expro.sample1){
  print('in check')
  #print(dim(expro.sample1))
  all=apply(expro.sample1, 2, function(x){
    return(as.numeric(as.character(x)))
  })
  row.names(all)=row.names(expro.sample1)
  colnames(all)=colnames(expro.sample1)
  nacount=apply(all, 2, function(x){
    return(sum(is.na(x)))
  })
  mincount=apply(all, 2, function(x){
    #return(min(table(x)))
    if(length(table(x))==1){
      return(0)
    }else if(length(table(x))==2&min(table(x))<3){
      return(0)
    }
    return(1)
  })
  a.inds=which(nacount<0.6*nrow(expro.sample1)&mincount>0)
  if(length(a.inds)==0){
    dat=data.frame(Random1=runif(nrow(expro.sample1),0,1),Random2=runif(nrow(expro.sample1),0,1))
  }else if(length(a.inds)==1){
    dat=data.frame(all[,a.inds],Random1=runif(nrow(expro.sample1),0,1))
  }else{
    dat=all[,a.inds]
  }
  row.names(dat)=row.names(expro.sample1)
  print('end check')
  return(dat)
}

tryCatch({
	  print('SUCC IN')
		outInfo=c(ncol(expro),nrow(expro),nclinical)
		
		expro=na.omit(expro)
		sds=apply(expro,1,sd)
		#print(cutSd)
		#print(length(which(sds>cutSd)))
		expro=expro[which(sds>cutSd),] ##########6/26修改cutsd########
		#print(head(expro))
		expro=t(expro)
		#print(head(expro))
		gsg = goodSamplesGenes(expro, verbose = 3); ##检测缺失�?#######
    #print(gsg$allOK)
		if (!gsg$allOK){
			# Optionally, print the gene and sample names that were removed:
			if (sum(!gsg$goodGenes)>0)
				printFlush(paste("Removing genes:",paste(names(expro)[!gsg$goodGenes], collapse = ",")));
			if (sum(!gsg$goodSamples)>0)
				printFlush(paste("Removing samples:",paste(rownames(dataExpr)[!gsg$goodSamples], collapse = ",")));
			# Remove the offending genes and samples from the data:
			expro = expro[gsg$goodSamples, gsg$goodGenes]
		}
		nGenes = ncol(expro)
		nSamples = nrow(expro)
		print(dim(expro))
		sampleTree = hclust(dist(expro), method = "average")
		if(cutHight==-1){
		  cutHight.tmp=3*sd(sampleTree$height)+mean(sampleTree$height)
		  tree.tmp=cutree(sampleTree,h = cutHight.tmp)
		  #which(tree.tmp==names(sort(-table(tree.tmp))[1]))
		  #print(table(tree.tmp))
		  #print(names(sort(table(tree.tmp),decreasing = T))[1])
		  N=sort(table(tree.tmp),decreasing = T)[1]
		  #print(cutHight.tmp)
		  #print(quantile(sampleTree$height,seq(0,1,0.05)))
		  if(N/nrow(expro)>0.8){
		    cutHight=cutHight.tmp
		  }else{
		    cutHight=quantile(sampleTree$height,seq(0,1,0.05))['100%']+0.001
		  }
		}

		pdf('SampleCluster.pdf',width = 8,height = 6)
		par(cex = 0.6)
		par(mar = c(0,4,2,0))
		plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,cex.axis = 1.5, cex.main = 2)
		if(cutHight!=0){
			abline(h = cutHight, col = "red");
		} 
		dev.off()
		jpeg('SampleCluster.jpg',width = 800,height = 600)
		par(cex = 0.6)
		par(mar = c(0,4,2,0))
		plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,cex.axis = 1.5, cex.main = 2)
		if(cutHight!=0){
			abline(h = cutHight, col = "red");
		}
		dev.off()
		
		#####add Sample#####
		expro.sample=OutSampleCluster(sampleTree,expro.sample)
		#################
		if(cutHight!=(-1)){
		  tree=cutree(sampleTree,h=cutHight)
		  data1=expro[which(tree==names(sort(table(tree),decreasing = T))[1]),]
		  #if(nclinical==1){
		    dat=data.frame(gb=expro.sample[which(tree==names(sort(table(tree),decreasing = T))[1]),])
		    colnames(dat)=colnames(expro.sample)
		    expro.sample=dat
		  #}else{
		  #  expro.sample=expro.sample[which(tree==order(-table(tree))[1]),]
		  #}
		}else{
			data1=expro
		}
		expro.sample=checkClinicalData(expro.sample)
		print(dim(data1))
		print(dim(expro.sample))
		write.table(cbind(row.names(data1),data1),file ="WGCNA_Input_data.txt",sep = '\t',quote = F,row.names = F)
		write.table(cbind(row.names(expro.sample),expro.sample),file ="WGCNA_Input_SampleInfo.txt",sep = '\t',quote = F,row.names = F)

		
		outInfo=c(outInfo,ncol(data1),nrow(data1))
		powers = c(c(1:10), seq(from = 12, to=20, by=2))
		sft = pickSoftThreshold(data1, powerVector = powers, verbose = 5,networkType = type)
		cutPower=sft$powerEstimate
		nSamples=nrow(data1)
		if (is.na(cutPower)){
		  cutPower = ifelse(nSamples<20, ifelse(type == "unsigned", 9, 18),
				ifelse(nSamples<30, ifelse(type == "unsigned", 8, 16),
				ifelse(nSamples<40, ifelse(type == "unsigned", 7, 14),
				ifelse(type == "unsigned", 6, 12)) 
					)
				)
		}
		scorePowers=-sign(sft$fitIndices[,3])*sft$fitIndices[,2]
		scoreIndices=sft$fitIndices[,5]
		cutScore=which(powers==cutPower)
		colors=rep('blue',length(powers))
		colors[cutScore]='red'

		pdf('PowerList.pdf',width = 8,height = 6)
		par(mfrow = c(1,2));
		plot(sft$fitIndices[,1], scorePowers,
			xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
			main = paste("Scale independence"));
		text(sft$fitIndices[,1], scorePowers,
		labels=powers,cex=cex1,col=colors);
		abline(h=0.85,col="red")
		# Mean connectivity as a function of the soft-thresholding power
		plot(sft$fitIndices[,1], sft$fitIndices[,5],
			xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
			main = paste("Mean connectivity"))
		text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col=colors)
		dev.off()
		jpeg('PowerList.jpg',width = 800,height = 600)
		par(mfrow = c(1,2));
		plot(sft$fitIndices[,1], scorePowers,
			xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
			main = paste("Scale independence"));
		text(sft$fitIndices[,1], scorePowers,
		labels=powers,cex=cex1,col=colors);
		abline(h=0.85,col="red")

		# Mean connectivity as a function of the soft-thresholding power
		plot(sft$fitIndices[,1], sft$fitIndices[,5],
			xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
			main = paste("Mean connectivity"))
		text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col=colors)
		dev.off()

		outInfo=c(outInfo,cutPower)

		#print("+==One-step network construction and module detection==")
		net = blockwiseModules(
			data1,
			power = cutPower,
			maxBlockSize = BlockSize,
			TOMType = type , minModuleSize = minModuleSize,
			reassignThreshold = 0, mergeCutHeight = mergeCutHeight,
			deepSplit = deepSplit,
			numericLabels = TRUE, pamRespectsDendro = FALSE,
			saveTOMs = T,
			corType = corType,
			saveTOMFileBase = paste0(ePath,".tom"),
			verbose = 3
		)
		print(table(net$colors))

		mergedColors = labels2colors(net$colors)
		write.table(cbind(colnames(data1),mergedColors)[order(mergedColors),],file = "GeneModuleClass.txt",sep = '\t',quote = F,row.names = F)
		write.table(table(mergedColors),file = "GeneModuleClassCount.txt",sep = '\t',quote = F,row.names = F)

		pdf('DendroAndColors.pdf',width = 8,height = 6)
		plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
			"Module colors",
			dendroLabels = FALSE, hang = 0.03,
			addGuide = TRUE, guideHang = 0.05)
		dev.off()
		#class(mergedColors[net$blockGenes[[1]]])
		jpeg('DendroAndColors.jpg',width = 800,height = 600)
		plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
			"Module colors",
			dendroLabels = FALSE, hang = 0.03,
			addGuide = TRUE, guideHang = 0.05)
		dev.off()
		print(length(table(mergedColors)))
		print(length(table(mergedColors))==1&&names(table(mergedColors))[1]=='grey')
    if(length(table(mergedColors))==1&&names(table(mergedColors))[1]=='grey'){
      jsTxt=rbind(
        paste0("var args={cutHight:",round(cutHight,2),",BlockSize:",BlockSize,",minModuleSize:",minModuleSize,",mergeCutHeight:",mergeCutHeight,",cutR:",cutR,",cutNet:",cutNet,",cutSd:",cutSd,",module:",0,"}")
        ,paste0("var info=[",paste0(outInfo,collapse = ","),"]")
        ,paste0("var mCount=[",paste0(table(mergedColors),collapse = ","),"]")
        ,paste0("var mName=['",paste0(names(table(mergedColors)),collapse = "','"),"']")
        ,paste0("var mHubCount=[]")
        ,paste0("var mHubName=[]")
        ,paste0("var mNetNodeCount=[]")
        ,paste0("var mNetEdgeCount=[]")
        ,paste0("var mNetName=[]")
      )
      write.table(jsTxt,file = "inputJs.js",quote = F,row.names = F,col.names = F)
      print('END')
    }
		else{
		
		MEs0 = moduleEigengenes(data1, mergedColors)$eigengenes
		MEs = orderMEs(MEs0);
		colnames(MEs)=gsub('^ME','',colnames(MEs))
		me.inds=which(colnames(MEs)!='grey')
		if(length(me.inds)==1){
		  cns=colnames(MEs)[me.inds]
		  MEs=as.data.frame(MEs[,me.inds])
		  colnames(MEs)=cns
		}else if(length(me.inds)>1){
		  MEs=MEs[,me.inds]
		}
		
		geneTraitSignificance = as.data.frame(cor(data1,expro.sample, use = "p"));
		GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nrow(MEs)));
		names(geneTraitSignificance) = paste("GS.", colnames(expro.sample), sep="");
		names(GSPvalue) = paste("p.GS.", colnames(expro.sample), sep="")
		geneModuleMembership = as.data.frame(cor(data1, MEs, use = "p"))
		MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
		names(geneModuleMembership) = paste("MM", names(geneModuleMembership), sep="");
		names(MMPvalue) = paste("p.MM", names(geneModuleMembership), sep="");
		write.table(cbind(Genes=row.names(geneTraitSignificance),geneTraitSignificance,GSPvalue,Module=mergedColors)
		    ,file = paste0("GeneTraitSignificance.txt"),sep = '\t',quote = F,row.names = F)
		sups=c('GeneTraitSignificance.txt')
		h=1
		w=1
		if(ncol(MEs)<=3){
		  w=ncol(MEs)
		}else{
		  mds=c(ncol(MEs)%%3,ncol(MEs)%%4,ncol(MEs)%%5)
		  sb=c(3,4,5)[which(mds==min(mds))[1]]
		  h=ceiling(ncol(MEs)/sb)
		  w=sb
		}
		for(gi in 1:ncol(geneTraitSignificance)){
		  sups=c(sups,paste0(colnames(expro.sample)[gi],'_ModuleMembership2GeneSignificance'))
		pdf(paste0(colnames(expro.sample)[gi],'_ModuleMembership2GeneSignificance.pdf'),width = w*3,height = h*3)
		par(mfrow=c(h,w))
		for(m in colnames(MEs)){
		module = m
		column = match(module, colnames(MEs));
		moduleGenes = mergedColors==module;
		verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
		                   abs(geneTraitSignificance[moduleGenes, gi]),
		                   xlab = paste("Module Membership in", module, "module"),
		                   ylab = "Gene significance",
		                   main = paste("Module membership vs. gene significance\n"),
		                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
		}
		dev.off()
		jpeg(paste0(colnames(expro.sample)[gi],'_ModuleMembership2GeneSignificance.jpg'),width = w*200,height = h*200)
		par(mfrow=c(h,w))
		for(m in colnames(MEs)){
		  module = m
		  column = match(module, colnames(MEs));
		  moduleGenes = mergedColors==module;
		  verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
		                     abs(geneTraitSignificance[moduleGenes, gi]),
		                     xlab = paste("Module Membership in", module, "module"),
		                     ylab = "Gene significance",
		                     main = paste("Module membership vs. gene significance\n"),
		                     cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
		}
		dev.off()
		}
		error.bar <- function(x, y, upper, lower=upper, length=0.1,...){
		       if(length(x) != length(y) | length(y) !=length(lower) | length(lower) != length(upper))
		         stop("vectors must be same length")
		       arrows(x,y+upper, x, y-lower, angle=90, code=3, length=length, ...)
		       }
		gsMeans=rbind()
		gsSds=rbind()
		for(m in colnames(MEs)){
		  tmp=abs(geneTraitSignificance[which(mergedColors==m),])
		  y.means <- apply(tmp,2,mean)
		  y.sd <- apply(tmp,2,sd)
		  gsSds=rbind(gsSds,y.sd)
		  gsMeans=rbind(gsMeans,y.means)
		}
		row.names(gsMeans)=colnames(MEs)
		print('plot DistributionGeneSignificance.pdf')
		#dev.off()
		pdf('DistributionGeneSignificance.pdf',width = ifelse(nrow(gsMeans)*0.6>6,nrow(gsMeans)*0.6,6)
		    ,height = ifelse(ncol(gsMeans)*2>6,ncol(gsMeans)*2,6))
		par(mfrow=c(ncol(gsMeans),1))
		print(gsMeans)
		for(gi in 1:ncol(gsMeans)){
		  mY=max(gsMeans[,gi]+1.96*as.numeric(gsSds[,gi])/10)
		  barx=barplot((gsMeans[,gi])
		               ,ylim = c(0,mY),ylab = 'Average gene significance'
		               ,xlab = paste0('Traits By ',colnames(expro.sample)[gi])
		               ,col = colnames(MEs))
		  error.bar(barx,as.numeric(gsMeans[,gi]), 1.96*as.numeric(gsSds[,gi])/10)
		}
		dev.off()
		print('end plot')
		jpeg('DistributionGeneSignificance.jpg',width = ifelse(nrow(gsMeans)*50>600,nrow(gsMeans)*50,600)
		     ,height = ifelse(ncol(gsMeans)*200>600,ncol(gsMeans)*200,600))
		par(mfrow=c(ncol(gsMeans),1))
		for(gi in 1:ncol(gsMeans)){
		  mY=max(gsMeans[,gi]+1.96*as.numeric(gsSds[,gi])/10)
		  barx=barplot((gsMeans[,gi])
		               ,ylim = c(0,mY),ylab = 'Average gene significance'
		               ,xlab = paste0('Traits By ',colnames(expro.sample)[gi])
		               ,col = colnames(MEs))
		  error.bar(barx,as.numeric(gsMeans[,gi]), 1.96*as.numeric(gsSds[,gi])/10)
		}
		dev.off()
		
		
		#dim(data1)
		#expro.sample
		moduleTraitCor = cor(MEs, expro.sample , use = "p")
		moduleTraitPvalue = corPvalueStudent(moduleTraitCor,nrow(MEs))
		textMatrix = paste(signif(moduleTraitCor, 2), "(",
			signif(moduleTraitPvalue, 1), ")", sep = "");
		dim(textMatrix) = dim(moduleTraitCor)
		pdf("ModuleTraitRelationships.pdf",width = 8,height = 6)
		par(mar = c(6, 8.5, 3, 3));
		# Display the correlation values within a heatmap plot
		labeledHeatmap(Matrix = moduleTraitCor,
			xLabels = names(expro.sample),
			yLabels = names(MEs),
			ySymbols = names(MEs),
			colorLabels = FALSE,
			colors = blueWhiteRed(255),
			textMatrix = textMatrix,
			setStdMargins = FALSE,
			cex.text = 0.5,
			zlim = c(-1,1),
			main = paste("Module-trait relationships"))
		dev.off()
		jpeg("ModuleTraitRelationships.jpg",width = 800,height = 600)
		# Display the correlation values within a heatmap plot
		labeledHeatmap(Matrix = moduleTraitCor,
			xLabels = names(expro.sample),
			yLabels = names(MEs),
			cex.lab = 0.5,
			ySymbols = names(MEs),
			colorLabels = FALSE,
			colors = blueWhiteRed(255),
			textMatrix = textMatrix,
			setStdMargins = FALSE,
            cex.text = 0.5,
            zlim = c(-1,1),
            main = paste("Module-trait relationships"))
		dev.off()
		#print(ncol(MEs))
		dim(moduleTraitPvalue)
		allTraitCor=rbind()
		for(i in 1:ncol(MEs)){
			#inds=which(moduleTraitPvalue[i,]<0.5&moduleTraitCor[i,]>0)
			#if(length(inds)>0){
		  inds=which(moduleTraitPvalue[i,]<1.1)
				mes=cbind(rep(colnames(MEs)[i],length(inds))
				          ,colnames(expro.sample)[inds]
				          ,moduleTraitCor[i,inds],moduleTraitPvalue[i,inds])
				allTraitCor=rbind(allTraitCor,mes)
			#}
		}
		if(!is.null(allTraitCor)&&nrow(allTraitCor)>0){
			colnames(allTraitCor)=c("Module","Trait","R","P.value")
			write.table(allTraitCor,file ="ModuleTraitRelationships.txt",sep = '\t',quote = F,row.names = F)
		}else{
			write.table(allTraitCor,file ="ModuleTraitRelationships.txt",sep = '\t',quote = F,row.names = F)
		}

		hubCount=rbind()
		for(m in colnames(MEs)){
			tmp=data1[,which(mergedColors==m)]
			moduleTraitCor = cor(MEs[,which(colnames(MEs)==m)], tmp , use = "p");
			moduleTraitPvalue = corPvalueStudent(moduleTraitCor,nrow(MEs))
			mgCor=data.frame(colnames(tmp),moduleTraitCor[1,],moduleTraitPvalue[1,])[order(-moduleTraitCor[1,]),]
			colnames(mgCor)=c("gene","R","P.value")
			write.table(mgCor,file = paste0("ModuleHubGeneDetail_",m,".txt"),sep = '\t',quote = F,row.names = F)
			#cutR=0.6
			inds=which(moduleTraitCor>cutR)
			if(length(inds)>0){
				hubCount=rbind(hubCount,c(m,length(inds)))    
				L=0
				W=0
				if(length(inds)<3){
					L=1
					W=length(inds)
				}else if(length(inds)==4){
					L=2
					W=2
				}else{
					m4=length(inds)%%4
					if(4-m4<m4) m4=4-m4
					m3=length(inds)%%3
					if(3-m3<m3) m3=3-m3
					if(m4<m3){
						L=floor((length(inds))/4)
						if(m4!=0) L=L+1
							W=4   
					}else{
						L=floor((length(inds))/3)
						if(m3!=0) L=L+1
							W=3   
					}
				}
    
				pdf(paste0("ModuleHubGene_",m,".pdf"),width = 3*W,height = 3*L)
				par(mfrow=c(L,W),mar = c(5, 5, 1, 1))
				for(i in inds){
					plot(MEs[,which(colnames(MEs)==m)],tmp[,i],xlab = paste0("Module Membersship in ",m," module")
						,ylab = paste0(colnames(tmp)[i]," Expression"),cex=0.5
						,ylim = c(0,quantile(tmp[,i],probs = seq(0, 1, 0.05))['95%'])
						,xlim = c(0,quantile(MEs[,which(colnames(MEs)==m)],probs = seq(0, 1, 0.05))['95%'])
						,cex.lab=0.8,col="blue"
					)
					label=paste("p=",round(moduleTraitPvalue[i],3))
					if(moduleTraitPvalue[i]<0.001) label="p<0.001"
					legend("bottomright", c(paste("R=",round(moduleTraitCor[i],2)), label))
				}
				dev.off()
			}
		}

		if(ncol(MEs)>2){
  		pdf('ModuleRelationships.pdf',width = 8,height = 8)
  		par(mar = c(60, 8.5, 3, 3))
  		plotEigengeneNetworks(MEs, ""
  		                      #, excludeGrey=F
  		                      ,marDendro = c(0, 6, 2, 4), marHeatmap = c(4,6, 1, 1),
  				cex.lab = 0.8, xLabelsAngle = 45)
  		dev.off()
  
  		jpeg('ModuleRelationships.jpg',width = 800,height = 800)
  		par(mar = c(60, 8.5, 3, 3))
  		plotEigengeneNetworks(MEs, "", marDendro = c(0, 6, 2, 4), marHeatmap = c(4,6, 1, 1),
  				cex.lab = 0.8, xLabelsAngle = 45)
  		dev.off()
    }
		TOM = TOMsimilarityFromExpr(data1, power = cutPower,corType=corType, networkType=type);

		netCount=rbind()
		for(m in unique(mergedColors)){
			#m=unique(mergedColors)[1]
			probes=colnames(data1)
			inModule=is.finite(match(mergedColors,m))
			modProbes=probes[inModule]
			modTOM = TOM[inModule, inModule]
			dim(modTOM)
			dimnames(modTOM) = list(modProbes, modProbes)
			# Export the network into edge and node list files for Cytoscape
			cyt = exportNetworkToCytoscape(modTOM,
					#edgeFile=paste("CytoEdge",paste(m,collapse="-"),".txt",sep=""),
                    #nodeFile=paste("CytoNode",paste(m,collapse="-"),".txt",sep=""),
                    weighted = TRUE, threshold = cutNet,nodeNames=modProbes,
                    nodeAttr = mergedColors[inModule])
			nodeD=table(c(as.character(cyt$edgeData[,1])
				,as.character(cyt$edgeData[,2])
			))
			edg1=data.frame(as.character(cyt$edgeData[,1]),as.character(cyt$edgeData[,2]),cyt$edgeData[,3])
			colnames(edg1)=c("Node1","Node2","Weight")
			nodes=data.frame(as.character(cyt$nodeData[,1])
				,as.character(cyt$nodeData[,3])
				,nodeD[match(as.character(cyt$nodeData[,1]),names(nodeD))])
			nodes=nodes[,-3]
			colnames(nodes)=c("Node","Module","Degree")
			nodes=nodes[order(-nodes[,3]),]  
			netCount=rbind(netCount,c(m,nrow(nodes),nrow(edg1)))
			write.table(edg1,file =paste("CytoEdge",m,".txt",sep=""),sep = '\t',quote = F,row.names = F)
			write.table(nodes,file =paste("CytoNode",m,".txt",sep=""),sep = '\t',quote = F,row.names = F)
		}
		colnames(netCount)=c("Module","NodeCount","EdgeCount")
		write.table(netCount,file ="Network_count.txt",sep = '\t',quote = F,row.names = F)

		jsTxt=rbind(
		paste0("var args={cutHight:",round(cutHight,2),",BlockSize:",BlockSize,",minModuleSize:",minModuleSize,",mergeCutHeight:",mergeCutHeight,",cutR:",cutR,",cutNet:",cutNet,",cutSd:",cutSd,",module:",ncol(MEs),"}")
			,paste0("var info=[",paste0(outInfo,collapse = ","),"]")
			,paste0("var mCount=[",paste0(table(mergedColors),collapse = ","),"]")
			,paste0("var mName=['",paste0(names(table(mergedColors)),collapse = "','"),"']")
  			,paste0("var mHubCount=[",paste0(hubCount[,2],collapse = ","),"]")
			,paste0("var mHubName=['",paste0(hubCount[,1],collapse = "','"),"']")
  			,paste0("var mNetNodeCount=[",paste0(netCount[,2],collapse = ","),"]")
			,paste0("var mNetEdgeCount=[",paste0(netCount[,3],collapse = ","),"]")
			,paste0("var mNetName=['",paste0(netCount[,1],collapse = "','"),"']")
		  ,paste0('var gs=["',paste0(colnames(expro.sample),collapse = '","'),'"]')
		)
		write.table(jsTxt,file = "inputJs.js",quote = F,row.names = F,col.names = F)
	}
},error=function(e){
	files=c()
	cat("ERROR :",conditionMessage(e),"\n")
})

if(length(files)>0){
	write.table(files,file = 'run.log.txt',sep = '\t',quote = F,row.names = F,col.names =F)
}


##LASSO Cox regression
surv <- read.csv('Patient_survivalData.csv')
sig_gene <- read.csv('Patient_expressionDatafile.csv', row.names = 1)
cvfit = cv.glmnet(t(sig_gene), Surv(surv$OS, surv$EVENT), 
                  family = "cox") 
coef.min = coef(cvfit, s = "lambda.min")
active.min = which(!coef.min == 0)
geneids <- rownames(sig_gene)[active.min]
index.min = coef.min[active.min]
new <- cbind(geneids, index.min) %>%
  data.frame
new$index.min <- as.character(new$index.min) %>%
  as.numeric()
signature <- as.matrix(t(sig_gene[geneids,])) %*% as.matrix(index.min) 
signature <- data.frame(signature)
names(signature) <- 'Risk_score'
write.csv(signature, 'Risk Score for each patient.csv')


##survival decision tree
library(rpart)
library(rpart.plot)
pfit <- rpart(Surv(OS, EVENT) ~ Age + Gender + p-stage + HRS, data = df)
rpart.plot(pfit)

