R/PWcTalkNW.R

In Haocan223/MetaGSCA: Meta Gene Set Co-expression Analysis

# PWcTalkNW(): to draw Pathway crosstalk diagram.
# INPUT: PWpair. *.PWpair.txt that was generated by PPIpairs.forPWpairs().
### Must include the three fields: "PW1" (name), "PW2" (name), and "minus.logP"
# INPUT: PWp. Two-column text file with no header. $metaP component of list output of dichotGSAR().
### 1st column: pathways; 2nd column: meta P-value
# OUTPUT: an igraph object "g", i.e., the PW crosstalking network.
# NOTE: Usually this function has to run at least two times. 1st time for an interactively tunable tkplot,
### and 2nd time for a PDF plot based on the layout of the tuned tkplot.
PWcTalkNW <- function(PWpair,PWp,layout=NULL,figname='PWcTalk',
  pdfW=10,pdfH=10,asp=0.7,vbase=15,ebase=2,vlbase=1,power=1/2)
{
  write.csv(PWpair,paste0('PWpair.',figname,'.csv'),row.names=FALSE,quote=FALSE)
  write.csv(PWp,paste0('PWp.',figname,'.csv'),row.names=FALSE,quote=FALSE)
	if (max(PWpair[,3])<=1)
		PWpair[,3] <- -log10(PWpair[,3])
	colnames(PWpair)[3] <- 'minus.logP'
	PW.w <- -log10(unlist(PWp[,2]))
  names(PW.w) <- PWp[,1]
	if (!require(igraph)) stop('igraph must be preinstalled!\n')
	g <- graph_from_data_frame(PWpair,directed=FALSE)
	V(g)$weight <- PW.w[V(g)$name]
	E(g)$weight <- PWpair$minus.logP
 	E(g)$color <- 'grey'
	if (is.null(layout)) {
    cat('Graph of PWcTalkNW has',vcount(g),'vertices and',ecount(g),'undirectional edges. See two CSV files saved on disk.\n')
  	cat('node weight (-log10(p)), min & max:\t',min(signif(PW.w,2)), '\t',max(signif(PW.w,2)),'\n')
  	cat('edge weight (-log10(p)), min & max:\t',min(signif(PWpair$minus.logP,2)),'\t',max(signif(PWpair$minus.logP,2)),'\n')
		tkid=tkplot(g,layout=layout_with_fr,vertex.color='lightgray',vertex.frame.color=NA,
      vertex.size=vbase*(V(g)$weight)^(power),vertex.label.cex=vlbase*degree(g)^(power),edge.width=ebase*(E(g)$weight)^(power))
 	  return(list(g=g,tkid=tkid))
	} else {
		pdf(paste(figname,'pdf',sep='.'),width=pdfW,height=pdfH)
    par(mar=par()$mar+c(0,10,0,10))
		plot.igraph(g,vertex.label=V(g)$name,asp=asp,layout=layout,
      vertex.size=vbase*(V(g)$weight)^(power),vertex.label.cex=vlbase*degree(g)^(power),
			vertex.label.color='black',vertex.color='lightgray',vertex.frame.color=NA,
			edge.color='lightgray',vertex.label.font=2,
			edge.width=ebase*(E(g)$weight)^(power))
		dev.off()
    return(list(g=g,tkid=NULL))
	}
}