# -------------
getNinSubCluster <- function(Cluster1){
	idx <- unique(Cluster1)
	n1 <- length(idx)
	vN <- matrix(1.0, n1,2)
	vN[,1] <- idx
	for (i in 1:n1){
		vN[i,2] <- sum(Cluster1 == idx[i])
	}
	return(vN)
}

# -------------
getTargetDrugBank <- function(iFileTmp){
	library(org.Hs.eg.db)
	#iFileTmp <- paste(dirDat0, '/uniprot links.csv', sep='')

	f1 <- iFileTmp
	drg.dbank <- read.delim(f1, header=T, sep=',');
	drg.dbank <- as.matrix(drg.dbank[,c(2,4)]);

	drg.dbank1 <- drg.dbank
	keysTmp <- unique(drg.dbank1[,2])
	#keysTmp <- drg.dbank1[,2]
	uniprotSym <- select(org.Hs.eg.db, keys=keysTmp, columns=c('UNIPROT', 'SYMBOL', 'ENTREZID'), keytype="UNIPROT")
	uniprotSym <- as.matrix(uniprotSym);
	uniprotSym1 <- uniprotSym[is.na(uniprotSym[,2]) != TRUE,]
	uniprotSym1 <- uniprotSym1[uniprotSym1[,2] != "",]

	drg.dbank2 <- matrix('test', 1,2);
	for (i in 1:dim(drg.dbank1)[1]){
		strt <- drg.dbank1[i,1];
		uidt <- drg.dbank1[i,2];
		symt <- uniprotSym1[uniprotSym1[,1] == uidt,2];
		if (length(symt)==0){print(i); next;}
		if (length(symt)==1){
			vt <- matrix('test', 1,2);
			vt[1] <- strt; vt[2]<-symt;
			drg.dbank2 <- rbind(drg.dbank2, vt);
		} else {
			for (j in 1:length(symt)){
				vt <- matrix('test', 1,2);
				vt[1] <- strt; vt[2]<-symt[j];
				drg.dbank2 <- rbind(drg.dbank2, vt);
			}
		}
	}

	drg.dbank2 <- drg.dbank2[2:dim(drg.dbank2)[1],]
	#f15 <- paste(dirDat0, 'allDrugTargetDrugBank.txt', sep='/');
	#write.table(drg.dbank2, f15, col.names=F, row.names=F, quote=F, sep='\t');
	return(drg.dbank2)
}


# -------------
getTargetStitch <- function(drg.target, drg.cid1){

	drg.cid <- drg.cid1[,2]
	drg.name <- drg.cid1[,1]
	drg.target <- drg.target[drg.target[,1] %in% drg.cid, ];
	drg.target1 <- drg.target[as.numeric(drg.target[,3]) > 700,]

	library(org.Hs.eg.db)
	x <- org.Hs.egSYMBOL; mapped_genes <- mappedkeys(x); xx <- as.list(x[mapped_genes]);  egID1 <- names(xx); symID <- xx; ### entrz id and gene symbol;
	xx <- as.list(org.Hs.egENSEMBLPROT2EG); ensgID <- names(xx); egID2 <- xx; ## entrz id and ensembl id;
	symID1 <- rep('test', length(egID1));
	for (i in 1:length(egID1)){
		symID1[i] <- symID[[i]];
	}

	kX1 <- 0;
	drg.target2 <- matrix('test', 1,4);
	idxt1 <- 1:length(ensgID);
	idxt2 <- 1:length(egID1);

	for (i in 1:dim(drg.target1)[1]){
		#print(i)
		cidt <- drg.target1[i,1];
		drg1 <- drg.name[which(drg.cid == cidt)]

		st <- drg.target1[i,3];
		ensgt <- drg.target1[i,2];
		ensgt <- substr(ensgt, 6,20);
		idt1 <- idxt1[ensgID %in% ensgt];

		if (length(idt1)<1){
			next
		}  # not found

		egidt <- egID2[[idt1]];
		if (length(egidt)==1){
			vt <- matrix('test', 1,4);
			idt2 <- idxt2[egID1 == egidt]
			vt[1]<-cidt; vt[2]<-symID1[idt2]; vt[3] <- st; vt[4] <- drg1
			if (length(drg1) != 1){print(i)}
			kX1 <- kX1 + 1
			drg.target2 <- rbind(drg.target2, vt);

		} else {
			if (length(egidt)>1){
				#print(i)
				#print(egidt)
				#print(length(egidt))  # have more than one entrezID corresonding to the ensembl id;
				for (j in 1:length(egidt)){
					vt <- matrix('test', 1,3);
					idt2 <- idxt2[egID1 == egidt[j]]
					vt[1]<-cidt; vt[2]<-symID1[idt2]; vt[3] <- st; vt[4] <- drg1;
					if (length(drg1) != 1){print(i)}
					kX1 <- kX1 + 1
					drg.target2 <- rbind(drg.target2, vt);
				}
			}
		}
	}

	return(drg.target2)
}

# -------------
getName2CID2 <- function(drg.name){  # convert pubchem compound ID (CID), e.g., from STITCH to compound names

	nDrg <- length(drg.name)
	library(webchem)
	drg.cid <- matrix('test', 1,2);
	drg.cid[1,1] <- c("drgName");
	drg.cid[1,2] <- c("drgCID");

	for (i in 1:nDrg){
		print(i)
		strt <- drg.name[i];
		idt <- get_cid(strt, verbose=F);
		idt <- idt[[1]];
		nt <- length(idt);
		if (nt >= 1){
			vt <- matrix('test', 1,2)
			vt[1,1] <- strt;
			vt[1,2] <- idt[1];
			drg.cid <- rbind(drg.cid, vt);

		}


		# if (nt > 1){
		# 	for (j in 1:nt){
		# 		t <- cid_compinfo(idt[j],verbose=F); t1 <- tolower(t$synonyms); t2 <- t1[t1 %in% strt];
		# 		if (length(t2)>0){
		# 			drg.cid[i] <- idt[j];
		# 			break;
		# 		}
		# 	}
		# }
	}

	nCid <- dim(drg.cid)[1]
	for (i in 2:nCid){
		vt <- drg.cid[i,2];
		if (is.na(vt)){
			next
		}
		nt <- nchar(vt);
		if (nt > 8){
			drg.cid[i,2] <- 'NA'
			next
		}
		xt <- 1*10^(8-nt);
		strt <- paste('CID', as.character(xt),vt, sep='');
		drg.cid[i,2] <- gsub('CID1', 'CID0', strt);
	}
	#drg.cid[66]<-'CID109568614'; drg.cid[106] <- 'CID105362420'; drg.cid[447] <- 'CID100072385';
	#rownames(drg.cid) <- drg.name;
	#f06 <- paste(dirDat0, 'cmapDrugCIDs.txt', sep='/');
	#write.table(drg.cid, f06, col.names=F, sep='\t', quote=F);
	return(drg.cid)
}


# -------------
getPatient2DrugReverseDis2 <- function(geneSetUp, geneSetDown, gSymZs, rankMatrix, n.probe.all, n.drg){ ###
	m.bigS <- matrix(-2.0, 1, n.drg)
	idxt <- 1:n.probe.all;
	probe.set.up <- idxt[gSymZs%in%geneSetUp];
	if (length(probe.set.up)>1 & sum(is.na(probe.set.up)==TRUE)<1){
		up.ks = apply(rankMatrix, 2, GSEA.EnrichmentScore2, gene.set=probe.set.up)
	} else {
		up.ks = 0.0;
	}
	probe.set.down <- idxt[gSymZs%in%geneSetDown];
	if (length(probe.set.down)>1 & sum(is.na(probe.set.down)==TRUE)<1){
		down.ks = apply(rankMatrix, 2, GSEA.EnrichmentScore2, gene.set=probe.set.down)
	} else {
		down.ks = 0.0;
	}
	m.bigS = down.ks - up.ks
	return(m.bigS)
}

# ---------------
getSynScoreV3 <- function(tar0, tar1, sigNet1, w1, w2){### targets of d0; targets of d1; disease signaling network;

	n0 <- length(tar0);
	n1 <- length(tar1);
	vSet0 <- V(sigNet1)$name

	## score 1: Target Centrality
	sc1 <- closeness(sigNet1, vSet0);  sc2 <- betweenness(sigNet1, vSet0, directed=F); sc3 <- page.rank(sigNet1, algo="prpack", vids=vSet0)$vector; ### for normalization
	sc01 <- closeness(sigNet1, tar0);  sc02 <- betweenness(sigNet1, tar0, directed=F); sc03 <- page.rank(sigNet1, algo="prpack", vids=tar0)$vector;
	sc11 <- closeness(sigNet1, tar1);  sc12 <- betweenness(sigNet1, tar1, directed=F); sc13 <- page.rank(sigNet1, algo="prpack", vids=tar1)$vector;
	sc01 <- (sc01-min(sc1))/(max(sc1)-min(sc1)); sc02 <- (sc02-min(sc2))/(max(sc2)-min(sc2)); sc03 <- (sc03-min(sc3))/(max(sc3)-min(sc3));          ###normalization
	sc11 <- (sc11-min(sc1))/(max(sc1)-min(sc1)); sc12 <- (sc12-min(sc2))/(max(sc2)-min(sc2)); sc13 <- (sc13-min(sc3))/(max(sc3)-min(sc3));
	scTar0 <- (sum(sc01)+sum(sc02)+sum(sc03))/3.0;
	scTar1 <- (sum(sc11)+sum(sc12)+sum(sc13))/3.0;

	sScore <- w1*scTar0 + w2*scTar1;
	return(sScore);
}


# -------------
getSynScore2 <- function(tar0, tar1, sigNet1){### targets of d0; targets of d1; disease signaling network;

	n0 <- length(tar0);
	n1 <- length(tar1);
	vSet0 <- V(sigNet1)$name

	## score 1: Target Centrality
	sc1 <- closeness(sigNet1, vSet0);  sc2 <- betweenness(sigNet1, vSet0, directed=F); sc3 <- page.rank(sigNet1, algo="prpack", vids=vSet0)$vector; ### for normalization
	sc01 <- closeness(sigNet1, tar0);  sc02 <- betweenness(sigNet1, tar0, directed=F); sc03 <- page.rank(sigNet1, algo="prpack", vids=tar0)$vector;
	sc11 <- closeness(sigNet1, tar1);  sc12 <- betweenness(sigNet1, tar1, directed=F); sc13 <- page.rank(sigNet1, algo="prpack", vids=tar1)$vector;
	sc01 <- (sc01-min(sc1))/(max(sc1)-min(sc1)); sc02 <- (sc02-min(sc2))/(max(sc2)-min(sc2)); sc03 <- (sc03-min(sc3))/(max(sc3)-min(sc3));          ###normalization
	sc11 <- (sc11-min(sc1))/(max(sc1)-min(sc1)); sc12 <- (sc12-min(sc2))/(max(sc2)-min(sc2)); sc13 <- (sc13-min(sc3))/(max(sc3)-min(sc3));
	scTar0 <- (sum(sc01)+sum(sc02)+sum(sc03))/3.0;
	scTar1 <- (sum(sc11)+sum(sc12)+sum(sc13))/3.0;

	sScore <- scTar0 + scTar1;
	return(sScore);
	return(scTar0);
	return(scTar1);

}

# ---------
plotDrugDuoNetV1 <- function(net0, tar1, tar2, tar0){  #plot the driver signaling network

	net1 <- unique(net0)
	nEdge1 <- dim(net1)[1]
	nEdge0 <- dim(net0)[1]
	w1 <- rep(0.0, nEdge1)

	for (i in 1:nEdge1){
		a1 <- net1[i,1]
		b1 <- net1[i,2]
		for (j in 1:nEdge0){
			a2 <- net0[j,1]
			b2 <- net0[j,2]
			if (a1 == a2 & b1 == b2){
				w1[i] <- w1[i]+1.0
			}
		}
	}

	library(igraph)
	g1 <- graph.edgelist(net1)
	col1 <- V(g1)$name  # node name


	tar <- tar1
	nTar <- length(tar)
	if (nTar > 0){
		for (i in 1:nTar){
			str1 <- tar[i]
			col1 <- gsub(str1, 'red', col1)
			# col1 <- gsub(str1, 'red', col1)
		}
	}

	tar <- tar2;
	nTar <- length(tar)
	if (nTar > 0){
		for (i in 1:nTar){
			str1 <- tar[i]
			col1 <- gsub(str1, 'green', col1)
			# col1 <- gsub(str1, 'red', col1)
		}
	}

	tar <- tar0;
	nTar <- length(tar)
	if (nTar > 0){
		for (i in 1:nTar){
			str1 <- tar[i]
			col1 <- gsub(str1, 'cyan', col1)
			# col1 <- gsub(str1, 'red', col1)
		}
	}

	idx1 <- which(col1 != 'green' & col1 != 'red' & col1 != 'cyan')
	col1[idx1] <- 'gray'  # background node color
	# print(col1)

	V(g1)$color <- col1
	V(g1)$frame.color <- 'white'
	# ref of plot in igraph: http://igraph.org/r/doc/plot.common.html
	plot(g1, vertex.label.cex = 0.75, vertex.size=3, edge.arrow.mode=0, edge.arrow.size = 0.2, vertex.frame.color = 'white', vertex.label.font = 10, edge.width = 0.2)

}

# -Please download (add the function here) and cite the code from: https://rdrr.io/bioc/EnrichmentBrowser/src/R/GSEA.1.0.R 
# GSEA.EnrichmentScore2 <- function(gene.list, gene.set, weighted.score.type = 0, correl.vector = NULL) {
#