# this function is a wrapper that computes and plots correlation matrices of cell types 
# it takes as input average cell type matrices and differentially expressed genes
# it calculates correlation significance with a permutation test (1000 permutations)
# requires the functions in the file Comparative_functions.R: 
# requires the package corrplot

# input:
# ExpressionTableSpecies1 = data.frame with the average expression data by cell type in Species 1 (e.g. output of the AverageExpression function in Seurat)
# DEgenesSpecies1 = vector with the names of the genes differentially expressed among the cell types of the Species 1 set
# ExpressionTableSpecies2 - data.frame with the average expression data by cell type in Species 2
# DEgenesSpecies2 = vector with the names of the genes differentially expressed among the cell types of the Species 2 set
# species1 = c('turtle','lizard','mouse')
# species2 = c('turtle','lizard','mouse','human','chicken')
# filename = to save plots (in pdf)
# Height and Width are plot parameters

# returns a plot and a list:
# - plot with correlation matrices calculated by taking (1) the intersection of DE genes in the two species, (2) the union of DE genes, (3) only the DE genes in Species 1, (4) only the DE genes in Species 2. 
# - plot with correlation matrices calculated by taking intersect, union, Species1 or Species2 differentially expressed transcription factors
#
# list with all the correlation matrices (intersect, intersect.TF, union, union.TF, species1, species1.TF, species2, species2.TF)
 

CorrComparePlot <- function(ExpressionTableSpecies1, DEgenesSpecies1, ExpressionTableSpecies2, DEgenesSpecies2, Species1, Species2, filename, Height=12, Width=10) {
	
	Method1 = "intersect"
	comp.intersect <- SpPermute(ExpressionTableSpecies1, species1=Species1, DEgenesSpecies1, ExpressionTableSpecies2, species2=Species2, DEgenesSpecies2, nPermutations=1000, genes.use= Method1,corr.method="spearman")
	comp_table.intersect <- t(comp.intersect[[1]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.intersect[[1]])])
	p_table.intersect <- t(comp.intersect[[3]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.intersect[[1]])])
	p_table.intersect <- (-p_table.intersect)
	
	 comp.intersect.TF <- SpPermuteTF(ExpressionTableSpecies1, species1=Species1, DEgenesSpecies1, ExpressionTableSpecies2, species2=Species2, DEgenesSpecies2, nPermutations=1000, genes.use= Method1,corr.method="spearman")
	 comp_table.intersect.TF <- t(comp.intersect.TF[[1]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.intersect.TF[[1]])])
	 p_table.intersect.TF <- t(comp.intersect.TF[[3]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.intersect.TF[[1]])])
	 p_table.intersect.TF <- (-p_table.intersect.TF)
	
	Method2 = "union"
	comp.union <- SpPermute(ExpressionTableSpecies1, species1=Species1, DEgenesSpecies1, ExpressionTableSpecies2, species2=Species2, DEgenesSpecies2, nPermutations=1000, genes.use= Method2,corr.method="spearman")
	comp_table.union <- t(comp.union[[1]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.union[[1]])])
	p_table.union <- t(comp.union[[3]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.union[[1]])])
	p_table.union <- (-p_table.union)
	
	 comp.union.TF <- SpPermuteTF(ExpressionTableSpecies1, species1=Species1, DEgenesSpecies1, ExpressionTableSpecies2, species2=Species2, DEgenesSpecies2, nPermutations=1000, genes.use= Method2,corr.method="spearman")
	 comp_table.union.TF <- t(comp.union.TF[[1]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.union.TF[[1]])])
	 p_table.union.TF <- t(comp.union.TF[[3]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.union.TF[[1]])])
	 p_table.union.TF <- (-p_table.union.TF)
	
	# Method3 = "species1"
	# comp.species1 <- SpPermute(ExpressionTableSpecies1, species1=Species1, DEgenesSpecies1, ExpressionTableSpecies2, species2=Species2, DEgenesSpecies2, nPermutations=1000, genes.use= Method3,corr.method="spearman")
	# comp_table.species1 <- t(comp.species1[[1]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.species1[[1]])])
	# p_table.species1 <- t(comp.species1[[3]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.species1[[1]])])
	# p_table.species1 <- (-p_table.species1)
	
	# comp.species1.TF <- SpPermuteTF(ExpressionTableSpecies1, species1=Species1, DEgenesSpecies1, ExpressionTableSpecies2, species2=Species2, DEgenesSpecies2, nPermutations=1000, genes.use= Method3,corr.method="spearman")
	# comp_table.species1.TF <- t(comp.species1.TF[[1]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.species1.TF[[1]])])
	# p_table.species1.TF <- t(comp.species1.TF[[3]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.species1.TF[[1]])])
	# p_table.species1.TF <- (-p_table.species1.TF)
	
	# Method4 = "species2"
	# comp.species2 <- SpPermute(ExpressionTableSpecies1, species1=Species1, DEgenesSpecies1, ExpressionTableSpecies2, species2=Species2, DEgenesSpecies2, nPermutations=1000, genes.use= Method4,corr.method="spearman")
	# comp_table.species2 <- t(comp.species2[[1]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.species1[[1]])])
	# p_table.species2 <- t(comp.species2[[3]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.species1[[1]])])
	# p_table.species2 <- (-p_table.species2)
	
	# comp.species2.TF <- SpPermuteTF(ExpressionTableSpecies1, species1=Species1, DEgenesSpecies1, ExpressionTableSpecies2, species2=Species2, DEgenesSpecies2, nPermutations=1000, genes.use= Method4,corr.method="spearman")
	# comp_table.species2.TF <- t(comp.species2.TF[[1]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.species2.TF[[1]])])
	# p_table.species2.TF <- t(comp.species2.TF[[3]][1:ncol(ExpressionTableSpecies1),(ncol(ExpressionTableSpecies1)+1):nrow(comp.species2.TF[[1]])])
	# p_table.species2.TF <- (-p_table.species2.TF)

pdf(paste0(filename, ".pdf"), height=Height, width = Width)

	col1 <- colorRampPalette(c("darkblue", "white","darkred"))
	par(mfrow=c(2,2), oma=c(0,0,2,0)+0.5)
	
	corrplot(comp_table.intersect, order="original",tl.pos="lt", method="color", tl.col="black",cl.lim=c(min(comp_table.intersect),max(comp_table.intersect)), is.corr=F,tl.cex=0.7, col=col1(200), sig.level=(-0.05),insig="pch", pch=19, pch.cex=0.25,pch.col="black",p.mat=p_table.intersect, main= paste(Method1,",",comp.intersect[[8]], "genes", sep=" "),mar=c(3,1,5,1),cl.align.text="l")
		
	corrplot(comp_table.union, order="original",tl.pos="lt", method="color", tl.col="black",cl.lim=c(min(comp_table.union),max(comp_table.union)), is.corr=F,tl.cex=0.7, col=col1(200), sig.level=(-0.05),insig="pch", pch=19, pch.cex=0.25,pch.col="black",p.mat=p_table.union, main= paste(Method2,",",comp.union[[8]], "genes", sep=" "),mar=c(3,1,5,1),cl.align.text="l")
	
#	corrplot(comp_table.species1, order="original",tl.pos="lt", method="color", tl.col="black",cl.lim=c(min(comp_table.species1),max(comp_table.species1)), is.corr=F,tl.cex=0.7, sig.level=(-0.05),insig="pch", pch=19, pch.cex=0.25,pch.col="black",p.mat=p_table.species1, col=col1(200), main= paste(Method3,",",comp.species1[[8]], "genes", sep=" "),mar=c(3,1,5,1),cl.align.text="l")  
	
#	corrplot(comp_table.species2, order="original",tl.pos="lt", method="color", tl.col="black",cl.lim=c(min(comp_table.species2),max(comp_table.species2)), is.corr=F,tl.cex=0.7, sig.level=(-0.05),insig="pch", pch=19, pch.cex=0.25,pch.col="black",p.mat=p_table.species2, col=col1(200), main= paste(Method4,",",comp.species2[[8]], "genes", sep=" "),mar=c(3,1,5,1),cl.align.text="l") 	
	
#	mtext(paste(length(DEgenesSpecies1), Species1, "genes", "and", length(DEgenesSpecies2),  Species2,"genes,","all genes", sep=" "),outer=T)
	
#	par(mfrow=c(2,2), oma=c(0,0,2,0) + 0.5)
	
	corrplot(comp_table.intersect.TF, order="original",tl.pos="lt", method="color", tl.col="black",cl.lim=c(min(comp_table.intersect.TF),max(comp_table.intersect.TF)), is.corr=F,tl.cex=0.7, sig.level=(-0.05),insig="pch", pch=19, pch.cex=0.25,pch.col="black",p.mat=p_table.intersect.TF, col=col1(200), main= paste(Method1,",",comp.intersect.TF[[8]], "genes", sep=" "),mar=c(3,1,5,1),cl.align.text="l") 
	
	corrplot(comp_table.union.TF, order="original",tl.pos="lt", method="color", tl.col="black",cl.lim=c(min(comp_table.union.TF),max(comp_table.union.TF)), is.corr=F,tl.cex=0.7, sig.level=(-0.05),insig="pch", pch=19, pch.cex=0.25,pch.col="black",p.mat=p_table.union.TF, col=col1(200), main= paste(Method2,",",comp.union.TF[[8]], "genes", sep=" "),mar=c(3,1,5,1),cl.align.text="l")
	
#	corrplot(comp_table.species1.TF, order="original",tl.pos="lt", method="color", tl.col="black",cl.lim=c(min(comp_table.species1.TF),max(comp_table.species1.TF)), is.corr=F,tl.cex=0.7, sig.level=(-0.05),insig="pch", pch=19, pch.cex=0.25,pch.col="black",p.mat=p_table.species1.TF, col=col1(200), main= paste(Method3,",",comp.species1.TF[[8]], "genes", sep=" "),mar=c(3,1,5,1),cl.align.text="l") 
	 
#	corrplot(comp_table.species2.TF, order="original",tl.pos="lt", method="color", tl.col="black",cl.lim=c(min(comp_table.species2.TF),max(comp_table.species2.TF)), is.corr=F,tl.cex=0.7, sig.level=(-0.05),insig="pch", pch=19, pch.cex=0.25,pch.col="black",p.mat=p_table.species2.TF, col=col1(200), main= paste(Method4,",",comp.species2.TF[[8]], "genes", sep=" "),mar=c(3,1,5,1),cl.align.text="l") 
		
#	mtext(paste(length(DEgenesSpecies1),  Species1, "genes","and", length(DEgenesSpecies2),  Species2,"genes,","transcription factors", sep=" "),outer=T)
	
dev.off()

# objects.to.return = list(comp.intersect, comp.intersect.TF, comp.union, comp.union.TF, comp.species1, comp.species1.TF, comp.species2, comp.species2.TF)
# names(objects.to.return)  = c("intersect","intersect.TF","union", "union.TF","species1","species1.TF","species2","species2.TF")
# return(objects.to.return)	
	
}
