R/pcc_calc.R

In oacar/pgsNetwork: This package is used for ref

#!/home/oma21/bin/Rscript --vanilla
library(igraph)
library(tidyverse)
#Rcpp::sourceCpp('src/calculate_pcc.cpp')
args=commandArgs(TRUE)
net_fname=args[1]
#net_fname='analysis/data/derived_data/pcc_calculations/intwedgelist.graphml'
#outfile='analysis/data/derived_data/pcc_calculations/rcpp_pcccalc22'#args[2]
outfile=args[2]
nodelistout=args[3]
net=read_graph(net_fname,'graphml')
#nn1=read_graph('')
a=as_adj(net,attr = 'weight')
#class(a)
a_nonw=as_adj(net)
a_nonw.dense <- as.matrix(a_nonw)
a.dense<-as.matrix(a)
a.dense.na <- a.dense
a.dense.na[a_nonw.dense==0] <- NA 
#pcc=calculate_pcc(a_nonw.dense,a.dense)
pcc_symmetric <- cor(a.dense.na,use='pairwise.complete.obs') 
#colnames(pcc_symmetric) <- V(net)$name
#rownames(pcc_symmetric) <- V(net)$name
write.table(pcc_symmetric,outfile,sep='\t',col.names=FALSE,row.names=F,quote=F)
write.table(V(net)$name,nodelistout,sep='\t',col.names=F,row.names=F,quote=F)

#corna=cor(a.dense.na,use='pairwise.complete.obs')
#"pairwiseCount"  <-
    #function (x, y=NULL,diagonal=TRUE){
#           if(is.null(y)) {n <- t(!is.na(x)) %*% (!is.na(x)) } else { n <- t(!is.na(x)) %*% (!is.na(y)) } 
#       if(!diagonal & is.null(y)) diag(n) <- NA
#              return(n) }