R/genetic_correlation.r

   1  #SNOPSIS
   2
   3  #runs genetic correlation analyis.
   4  #correlation coeffiecients are stored in tabular and json formats
   5
   6  #AUTHOR
   7  # Isaak Y Tecle (iyt2@cornell.edu)
   8
   9
  10 options(echo = FALSE)
  11
  12 library(gplots)
  13 library(ltm)
  14 library(plyr)
  15 library(rjson)
  16
  17
  18 allargs<-commandArgs()
  19
  20 geneticDataFile <- grep("combined_gebvs",
  21                         allargs,
  22                         ignore.case=TRUE,
  23                         perl=TRUE,
  24                         value=TRUE
  25                       )
  26
  27 correTableFile <- grep("genetic_corre_table",
  28                        allargs,
  29                        ignore.case=TRUE,
  30                        perl=TRUE,
  31                        value=TRUE
  32                        )
  33
  34 correJsonFile <- grep("genetic_corre_json",
  35                       allargs,
  36                       ignore.case=TRUE,
  37                       perl=TRUE,
  38                       value=TRUE
  39                       )
  40
  41 geneticData <- read.table(geneticDataFile,
  42                           header = TRUE,
  43                           row.names = 1,
  44                           sep = "\t",
  45                           na.strings = c("NA"),
  46                           dec = "."
  47                           )
  48
  49 coefpvalues <- rcor.test(geneticData,
  50                          method="pearson",
  51                          use="pairwise"
  52                          )
  53
  54
  55 coefficients <- coefpvalues$cor.mat
  56 allcordata   <- coefpvalues$cor.mat
  57
  58 allcordata[lower.tri(allcordata)] <- coefpvalues$p.values[, 3]
  59 diag(allcordata) <- 1.00
  60
  61 pvalues <- as.matrix(allcordata)
  62 pvalues <- round(pvalues,
  63                  digits=2
  64                  )
  65
  66 coefficients <- round(coefficients,
  67                       digits=3
  68                       )
  69
  70 allcordata <- round(allcordata,
  71                     digits=3
  72                     )
  73
  74 #remove rows and columns that are all "NA"
  75 if ( apply(coefficients,
  76            1,
  77            function(x)any(is.na(x))
  78            )
  79     ||
  80     apply(coefficients,
  81           2,
  82           function(x)any(is.na(x))
  83           )
  84     )
  85   {
  86
  87     coefficients<-coefficients[-which(apply(coefficients,
  88                                             1,
  89                                             function(x)all(is.na(x)))
  90                                       ),
  91                                -which(apply(coefficients,
  92                                             2,
  93                                             function(x)all(is.na(x)))
  94                                       )
  95                                ]
  96   }
  97
  98
  99 pvalues[upper.tri(pvalues)]           <- NA
 100 coefficients[upper.tri(coefficients)] <- NA
 101
 102 coefficients2json <- function(mat){
 103     mat <- as.list(as.data.frame(t(mat)))
 104     names(mat) <- NULL
 105     toJSON(mat)
 106 }
 107
 108 traits <- colnames(coefficients)
 109
 110 correlationList <- list(
 111                    "traits"=toJSON(traits),
 112                    "coefficients"=coefficients2json(coefficients)
 113                    )
 114
 115 correlationJson <- paste("{",paste("\"", names(correlationList), "\":", correlationList, collapse=","), "}")
 116
 117 write.table(coefficients,
 118       file=correTableFile,
 119       col.names=TRUE,
 120       row.names=TRUE,
 121       quote=FALSE,
 122       dec="."
 123       )
 124
 125 write.table(correlationJson,
 126       file=correJsonFile,
 127       col.names=FALSE,
 128       row.names=FALSE,
 129       )
 130
 131 q(save = "no", runLast = FALSE)