cgi-bin/phenome/correlation.r

   1  #SNOPSIS
   2
   3  #commands for running correlation analyis,
   4  #and generating heatmap and the correlation
   5  #coefficients and their p-values
   6
   7
   8  #AUTHOR
   9  # Isaak Y Tecle (iyt2@cornell.edu)
  10
  11
  12 options(echo = FALSE)
  13
  14 library(gplots)
  15 library(RColorBrewer)
  16 library(ltm)
  17
  18 allargs<-commandArgs()
  19
  20 phenodata<-grep("phenodata",
  21                allargs,
  22                ignore.case=TRUE,
  23                perl=TRUE,
  24                value=TRUE
  25                )
  26
  27 cortable<-grep("corre_table",
  28                allargs,
  29                ignore.case=TRUE,
  30                perl=TRUE,
  31                value=TRUE
  32                )
  33
  34 print(c("correlation table file:", cortable))
  35
  36 heatmap<-grep("heatmap",
  37                allargs,
  38                ignore.case=TRUE,
  39                perl=TRUE,
  40                value=TRUE
  41                )
  42
  43 print(c("heatmap file:", heatmap))
  44
  45 #reading phenotype data into an R object
  46 phenodata<-read.csv(phenodata,
  47                     header=TRUE,
  48                     dec=".",
  49                     sep=",",
  50                     na.strings=c("NA", "-")
  51                    )
  52
  53 phenodata$ID=NULL
  54
  55 #running Pearson correlation analysis
  56 coefpvalues<-rcor.test(phenodata,
  57                        method="pearson",
  58                        use="pairwise"
  59                       )
  60
  61 coefficients<-coefpvalues$cor.mat
  62 allcordata<-coefpvalues$cor.mat
  63 allcordata[lower.tri(allcordata)]<-coefpvalues$p.values[, 3]
  64 diag(allcordata)<-1.00
  65
  66
  67 pvalues<-as.matrix(allcordata)
  68
  69
  70 pvalues<-round(pvalues,
  71                digits=2
  72                )
  73
  74
  75 #rounding correlation coeficients into 2 decimal places
  76 coefficients<-round(coefficients,
  77                     digits=2
  78                    )
  79
  80 allcordata<-round(allcordata,
  81                   digits=2
  82                   )
  83
  84 #remove rows and columns that are all "NA"
  85 if ( apply(coefficients,
  86            1,
  87            function(x)any(is.na(x))
  88            )
  89     ||
  90     apply(coefficients,
  91           2,
  92           function(x)any(is.na(x))
  93           )
  94     )
  95   {
  96
  97     coefficients<-coefficients[-which(apply(coefficients,
  98                                             1,
  99                                             function(x)all(is.na(x)))
 100                                       ),
 101                                -which(apply(coefficients,
 102                                             2,
 103                                             function(x)all(is.na(x)))
 104                                       )
 105                                ]
 106   }
 107
 108 print("correlation coeffients:")
 109 print(coefficients)
 110
 111 pvalues[upper.tri(pvalues)]<-NA
 112 coefficients[upper.tri(coefficients)]<-NA
 113
 114 png(file=heatmap,
 115     height=600,
 116     width=600,
 117     bg="transparent"
 118     )
 119
 120
 121 heatmap.2(coefficients,
 122           Rowv=FALSE,
 123           Colv=FALSE,
 124           dendrogram="none",
 125           na.rm=TRUE,
 126           col =  rev(colorRampPalette(brewer.pal(10,"Spectral"))(128)),
 127           trace="none",
 128           ColSideColors,
 129           RowSideColors,
 130           keysize = 1,
 131           density.info="none",
 132           lmat=rbind( c(0,3), c(2,1), c(0,4)),
 133           lhei=c(0.25, 4, 0.75),
 134           lwid=c(0.25, 4),
 135           cexRow = 1.25,
 136           cexCol = 1.25,
 137           margins = c(10, 6)
 138           )
 139
 140 dev.off()
 141
 142 write.table(coefficients,
 143       file=cortable,
 144       col.names=TRUE,
 145       row.names=TRUE,
 146       quote=FALSE,
 147       dec="."
 148       )
 149
 150 q(save = "no", runLast = FALSE)