Network analysis-next step after calculating correlation in R
1
0
Entering edit mode
7.4 years ago
mahnazkiani ▴ 50

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, schmid.m@access.uzh.ch 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

RNA-Seq R • 3.3k views
ADD COMMENT
1
Entering edit mode
7.4 years ago

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.

ADD COMMENT

Login before adding your answer.

Traffic: 1911 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6