gonad.R

# content of main.R
source("CSI.R")

#Cenix
data <- read.delim("../../Cenix_Digimap_raw_040820.csv", header = TRUE, row.names = 1, sep = "\t")
input_m <- data[,2:46]
phenoSig  <- as.matrix(input_m) 
phenoSig[ which(phenoSig <=1) ] = 0 
phenoSig[ which(phenoSig >1) ] = 1 

names <- row.names(phenoSig)
for (i in 1:length(names)){
	if ( substr(names[i], nchar(names[i]), nchar(names[i]) ) == "a" || substr(names[i], nchar(names[i]), nchar(names[i]) ) == "b" ){
		names[i] <- substr(names[i], 1, nchar(names[i])-1)
	}
}
row.names(phenoSig) <- names

sumPhenoSig <- rowSums(phenoSig)
filteredPhenoSig <- phenoSig[sumPhenoSig>0,]
tphenoSig <- t(filteredPhenoSig)
PCC <- cor(tphenoSig, method="pearson")
pccCSI <- CSI(PCC, 0.05)

if (exists ("result") ){
	rm(result)
}

for (i in 1:(nrow(pccCSI)-1) ){
	for (j in (i+1):nrow(pccCSI) ){
		if (!exists ("result") ){
			result <- c(row.names(filteredPhenoSig)[i], row.names(filteredPhenoSig)[j], pccCSI[i,j] )
		}else{
			result <- rbind(result, c(row.names(filteredPhenoSig)[i], row.names(filteredPhenoSig)[j], pccCSI[i,j] ) )
		}	
	}
}

write.table(result, file="../results/Cenix_2_CSI.txt", sep="\t", row.names = FALSE, quote =FALSE, col.names=FALSE)


#Gonad
phenoSig <- read.table("../../Greenetal_geneBinary2.txt", sep="\t", header=T, row.names=1)
phenoSig  <- as.matrix(phenoSig)
sumPhenoSig <- rowSums(phenoSig)
filteredPhenoSig <- phenoSig[sumPhenoSig>0,]
tphenoSig <- t(filteredPhenoSig)
PCC <- cor(tphenoSig, method="pearson")
pccCSI <- CSI(PCC, 0.05)

if (exists ("result") ){
	rm(result)
}
for (i in 1:(nrow(pccCSI)-1) ){
	for (j in (i+1):nrow(pccCSI) ){
		if (!exists ("result") ){
			result <- c(row.names(filteredPhenoSig)[i], row.names(filteredPhenoSig)[j], pccCSI[i,j] )
		}else{
			result <- rbind(result, c(row.names(filteredPhenoSig)[i], row.names(filteredPhenoSig)[j], pccCSI[i,j] ) )
		}	
	}
}

write.table(result, file="../results/gonad.txt", sep="\t", row.names = FALSE, quote =FALSE, col.names=FALSE)

# content of CSI.R

CSI <- function(input_m, offset){
#PCC (input_m) + network filter
#offset = 0.05 in the Green et al. paper

	geneNum <- length(input_m[,1])

	out_m <- matrix(data = 0, ncol=geneNum, nrow = geneNum, dimnames = list(row.names(input_m),colnames(input_m)))
	
	#transfer to CSI
	for (i in 1:geneNum){
		for (j in i:geneNum){
			cUvalueMinusCutOff = input_m[i,j] - offset
			x <- which(input_m[i,] >= cUvalueMinusCutOff )
			y <- which(input_m[j,] >= cUvalueMinusCutOff )
			z <- unique(c(x,y))
			####This calculation remove i,j. This will make the maximum of CSI = 1
			#out_m[i,j] <- 1- ( (length(z)-2) / geneNum)  #calculate specificity
			####This calculation do not remove i,j. Therefore, the maximum CSI will equals to 1-(2/geneNum). According to the paper figure, 1 - (# Genes Connected to A or B with PCC ≥ PCCAB- offset / Total # Genes in Screen). Therefore, I made this as the official version in the calculation.		
			out_m[i,j] <- 1- (length(z) / geneNum)  #calculate specificity
			out_m[j,i] <- out_m[i,j]
		}
	}
	return (out_m)

}