Network analysis-next step after calculating correlation in R

Question

[mahnazkiani]

I went through the R script that was online from the

http://www.pb.ethz.ch/downloads/

###############################################################################################################
# t_GCN.r script for computing targeted gene co-expression networks                                           #
#     1. "guide-query GCN" --- computes correlations between guide and query genes                            #
#         We are grateful to Marc W. Schmid, UZH, [e-mail] for improving the speed of this step #
#                                                                                                             #
#     2. "guide-query groups GCN" --- computes correlations between guide genes and groups of query genes     #
###############################################################################################################

# Required packages and libraries
source("http://bioconductor.org/biocLite.R")
library(Biobase)
require(multtest)
require(xtable)

################################################################################
# 1. "guide-query GCN" --- computes correlations between guide and query genes #
################################################################################


# set the threshold for P-values
alpha <- 0.05

# read in the file containing the expression data for the guide genes
guide_expression <-read.csv(file="T:/Your_WorkingDirectory/guide_expression.csv",header=TRUE, row.names=1)

guide_expression_names <- colnames(guide_expression)

# read in the file containing the expression data for the query genes
query_expression <- read.csv(file="T:/Your_WorkingDirectory/query_expression.csv",header=TRUE, row.names=1)

m <- ncol(query_expression)
n <- nrow(query_expression)

cornames <- colnames(query_expression)

data <- cbind(guide_expression, query_expression)

# calculate the Pearson correlation coefficients between guide genes and each group of query genes
# to use other similarity measures change for example "pearson" to "spearman"
corData <- cor(data, method ="pearson")[1:ncol(guide_expression),(ncol(guide_expression)+1):(m+ncol(guide_expression))]

# asses the statistical significance of positive Pearson correlation coefficients
all.pValue <- apply(corData, 1, function(x) 1-pnorm(sqrt(n-3)*0.5*(log(1+x)-log(1-x))))

# to assess the statistical significance of negative Pearson correlation coefficients run the entire script with the
# following modificantion: all.pValue <-apply(corData, 1, function(x) pnorm(sqrt(n-3)*0.5*(log(1+x)-log(1-x))))


all.pValue <- as.vector(all.pValue)

# Apply multiple testing correction (Bonferroni, Holm, FDR)

  bonf <- mt.rawp2adjp(all.pValue, proc="Bonferroni")
  pValue.bonf <- bonf$adjp[,2][order(bonf$index)]
  pValue.bonf <- matrix(pValue.bonf,nrow=m)
  pValue.bonf <-rbind(pValue.bonf,colSums((pValue.bonf<=alpha)))

  holm <- mt.rawp2adjp(all.pValue, proc="Holm")
  pValue.holm <- holm$adjp[,2][order(holm$index)]
  pValue.holm <- matrix(pValue.holm,nrow=m)
  pValue.holm <-rbind(pValue.holm,colSums((pValue.holm<=alpha)))

  fdr <- mt.rawp2adjp(all.pValue, proc="BH")
  pValue.fdr <- fdr$adjp[,2][order(fdr$index)]
  pValue.fdr <- matrix(pValue.fdr,nrow=m)
  pValue.fdr <-rbind(pValue.fdr,colSums((pValue.fdr<=alpha)))

  colnames(pValue.bonf) <- colnames(pValue.holm) <- colnames(pValue.fdr) <- guide_expression_names
  rownames(pValue.bonf) <- rownames(pValue.holm) <- rownames(pValue.fdr) <- c(cornames,"nr. signif.")

# save the results as HTML tables in your working directory)

  pValue.bonf <- xtable(pValue.bonf,digits=4)
    print(pValue.bonf, type="html", file=paste('guide_query_GCN',"bonf",".html",sep=""))

  pValue.holm <- xtable(pValue.holm,digits=4)
    print(pValue.holm, type="html", file=paste('guide_query_GCN',"holm",".html",sep=""))

  pValue.fdr <- xtable(pValue.fdr,digits=4)
    print(pValue.fdr, type="html", file=paste('guide_query_GCN',"fdr",".html",sep=""))

Now I have a .html file with correlation, how I can use this file to create a network. I tried to import it to Cytoscape but It didn't work. appreciate any help.

Thanks,
Mahnaz

Answer 1

[Jean-Karim Heriche]

As far as I know, you can't import html directly into cytoscape. It looks like you've got a correlation matrix. You can think of it as the adjacency matrix of a weighted undirected graph. Export it in a format that cytoscape can read e.g. sif. Depending on what you want to do/visualize, you can also filter the edges based on the associated p-values.