Running CyberT from R # Purpose You've got some microarray data that you would like to analyze with [CyberT](http://cybert.ics.uci.edu). The CyberT website provides an interface for running the analysis, but it also allows the option of downloading the [R](http://cran.r-project.org) code needed to perform the analysis on your own computer. R is an cross-platform, open-source, powerful statistics program which is mainly interacted with via a command-line. This is in contrast to other statistics software, such as SAS's JMP, which provide more menu-driven graphical interfaces for analysis. One could debate the relative merits of these approaches, but I will almost always provide tutorials for R because: it is available to everyone and includes more state-of-the-art techniques (such as CyberT analysis). See one my previous post on [Teaching Yourself R](http://anotherbio.tumblr.com/post/1095464070/teaching-yourself-r) for some basic tutorials. The purpose of this tutorial is show how to run CyberT in R on your own computer. # Pre-Requisites * An installation of R. Go to the [CRAN site](http://cran.r-project.org) for details. * The CyberT source files. The source code is available [here](http://cybert.ics.uci.edu/cgi-bin/CyberTReg-8.0.form.pl). It is licensed for non-commercial academic use. There are several sources available. I will be discussing the 'bayesreg' version of the code. * Some data to analyze as a text file in some delimited format. For purposes of the tutorial, let's use the example [CyberT Control + Experimental data](http://cybert.microarray.ics.uci.edu//sampleData/CyberT_C+E_DataSet) (which is exactly the data from the top of [this page](http://cybert.ics.uci.edu/cgi-bin/CyberT-8.0.form.pl?DATATYPE=CE)). This is a comma-delimited data set. # Tutorial ## Setup First, set up a new folder with the source file and your data file. I will refer to the new folder by the path /Users/mkayala/develop/CyberT/ and we'll name the source file 'bayesreg.R' and the data file 'testData.csv'. Then make a script file to save your R commands. R is an interpreted language with a console interface. We could work directly within the console, but then if you want to run or modify the analysis at a later date, you'll be stuck having to reenter everything from scratch. I use either [Emacs](http://www.gnu.org/software/emacs/) with [ESS](http://ess.r-project.org/) where you can send one line at a time to the R interpreter, or TextMate with applescript (subject of a future post) to do something similar. If you are unfamiliar with either of these, then I suggest copying and pasting code from the analysis code file to the R interpreter. Finally fire up your R interpreter. When I talk about a line of R code, I mean write it in the analysis file then send it to the R interpreter. ## Line by line First set the working directory of the R session. setwd('/Users/mkayala/develop/CyberT/') Then lets source the bayesreg.R file: source('bayesreg.R') And read in the data: data <- read.csv('testData.csv', header=TRUE) Then take a look at the data. I like to look at the first couple of lines, the dimensions etc: dim(data) head(data) tail(data) This should yield: > dim(data) [1] 2758 9 > head(data) X..gene.name C1 C2 C3 C4 E1 E2 E3 E4 1 gdhA 8.16e-05 7.59e-05 1.07e-04 1.02e-04 2.82e-04 2.45e-04 0.000296 2.70e-04 2 uvrA 1.27e-03 1.29e-03 1.27e-03 1.29e-03 1.03e-03 1.05e-03 0.001080 1.00e-03 3 hdhA 1.73e-06 1.39e-05 6.63e-06 2.96e-05 1.96e-04 2.32e-04 0.000185 2.42e-04 4 yecI 3.29e-05 3.21e-05 2.61e-05 3.32e-05 7.64e-05 8.07e-05 0.000092 9.02e-05 5 galP 3.64e-04 3.49e-04 4.00e-04 3.85e-04 2.10e-04 2.28e-04 0.000228 1.77e-04 6 malE 3.39e-04 3.34e-04 3.82e-04 3.69e-04 1.70e-04 2.14e-04 0.000218 2.00e-04 > tail(data) X..gene.name C1 C2 C3 C4 E1 E2 E3 E4 2753 b0817 5.52e-07 1.96e-05 2.86e-05 3.66e-05 2.61e-05 4.40e-05 4.02e-05 3.06e-05 2754 b2084 4.55e-05 4.58e-05 1.20e-06 1.66e-05 3.08e-05 4.43e-05 2.83e-05 3.76e-05 2755 citF 1.15e-04 8.48e-05 1.54e-04 1.27e-04 8.14e-05 7.79e-05 7.89e-05 8.26e-05 2756 yohL 3.78e-05 5.83e-05 1.39e-05 2.05e-05 6.97e-05 8.15e-05 7.39e-05 7.55e-05 2757 b0669 2.32e-04 1.86e-04 8.26e-06 5.72e-05 2.30e-04 2.28e-04 2.32e-04 1.76e-04 2758 yegH 7.03e-04 6.59e-04 5.51e-04 5.58e-04 6.38e-04 6.27e-04 6.42e-04 6.40e-04 Here we see that data is a data.frame with 2758 rows and 9 columns. The first column is really a name for the rows. Let's set the rownames attribute of the data frame, and remove this column from the actual 'data'. rownames(data) <- data[, 1] data <- data[, -1] Then take a look at the data again: head(data) > head(data) C1 C2 C3 C4 E1 E2 E3 E4 gdhA 8.16e-05 7.59e-05 1.07e-04 1.02e-04 2.82e-04 2.45e-04 0.000296 2.70e-04 uvrA 1.27e-03 1.29e-03 1.27e-03 1.29e-03 1.03e-03 1.05e-03 0.001080 1.00e-03 hdhA 1.73e-06 1.39e-05 6.63e-06 2.96e-05 1.96e-04 2.32e-04 0.000185 2.42e-04 yecI 3.29e-05 3.21e-05 2.61e-05 3.32e-05 7.64e-05 8.07e-05 0.000092 9.02e-05 galP 3.64e-04 3.49e-04 4.00e-04 3.85e-04 2.10e-04 2.28e-04 0.000228 1.77e-04 malE 3.39e-04 3.34e-04 3.82e-04 3.69e-04 1.70e-04 2.14e-04 0.000218 2.00e-04 Now we want to make a call to the 'bayesT' function to run the CyberT analysis. To see the first few lines of code that make up an R function, we can call head on it: head(bayesT) > head(bayesT) 1 function (aData, numC, numE, ppde = TRUE, betaFit = 1, bayes = TRUE, 2 winSize = 101, conf = 10) 3 { 4 if ((ceiling((winSize - 1)/2)) != ((winSize - 1)/2)) 5 stop("ERROR: winSize must be an odd number.") 6 numGene <- nrow(aData) All this does is show me what the arguments to the function are. These all correspond to the different form components on the web interface. The first is the data frame to work with. The second and third are the number of control and experimental columns, respectively. The fourth is whether or not to run a [PPDE](http://cybert.microarray.ics.uci.edu//help/index.html#ppde) analysis, and the fifth is the number of beta distributions to fit in the PPDE. Finally the last three have to do with the Bayesian part of the t-tests, whether or not to run, what window size to look around each gene, and the confidence in the windowed prior respectively. We'll stick with the defaults for all of these. However, play around with the different options to see what happens on our own. results <- bayesT(data, 4, 4); This should run and print out a few things about the PPDE fit. Then you have the results. head(results) > head(results) C1 C2 C3 C4 E1 E2 E3 E4 nC nE gdhA 8.16e-05 7.59e-05 1.07e-04 1.02e-04 2.82e-04 2.45e-04 0.000296 2.70e-04 4 4 uvrA 1.27e-03 1.29e-03 1.27e-03 1.29e-03 1.03e-03 1.05e-03 0.001080 1.00e-03 4 4 hdhA 1.73e-06 1.39e-05 6.63e-06 2.96e-05 1.96e-04 2.32e-04 0.000185 2.42e-04 4 4 yecI 3.29e-05 3.21e-05 2.61e-05 3.32e-05 7.64e-05 8.07e-05 0.000092 9.02e-05 4 4 galP 3.64e-04 3.49e-04 4.00e-04 3.85e-04 2.10e-04 2.28e-04 0.000228 1.77e-04 4 4 malE 3.39e-04 3.34e-04 3.82e-04 3.69e-04 1.70e-04 2.14e-04 0.000218 2.00e-04 4 4 meanC meanE stdC stdE fold rasdC rasdE gdhA 9.1625e-05 2.7325e-04 1.518560e-05 2.162368e-05 2.982265 3.768212e-05 5.769815e-05 uvrA 1.2800e-03 1.0400e-03 1.154701e-05 3.366502e-05 -1.230769 3.328681e-04 1.704221e-04 hdhA 1.2965e-05 2.1375e-04 1.216491e-05 2.752423e-05 16.486695 1.205830e-05 4.756146e-05 yecI 3.1075e-05 8.4825e-05 3.349005e-06 7.491495e-06 2.729686 1.705217e-05 3.274989e-05 galP 3.7450e-04 2.1075e-04 2.251666e-05 2.404683e-05 -1.776987 9.652963e-05 4.855883e-05 malE 3.5600e-04 2.0050e-04 2.322355e-05 2.174856e-05 -1.775561 8.488804e-05 4.590555e-05 bayesSDC bayesSDE bayesT bayesDF varRatio pVal ppde<(p) ppde(p) gdhA 3.522692e-05 5.376920e-05 5.650959 24 2.329795 8.085361e-06 0.9990119 0.9977560 uvrA 3.039204e-04 1.564813e-04 -1.404169 24 3.772200 1.730778e-01 0.7704639 0.5175552 hdhA 1.257637e-05 4.554638e-05 8.498692 24 13.115873 1.064495e-08 0.9999759 0.9999453 yecI 1.565623e-05 3.013016e-05 3.165952 24 3.703636 4.170175e-03 0.9682261 0.9301230 galP 8.883538e-05 4.592961e-05 -3.274795 24 3.740989 3.202789e-03 0.9724916 0.9392718 malE 7.835698e-05 4.329376e-05 -3.474011 24 3.275705 1.964593e-03 0.9789597 0.9532642 ROC(x) ROC(y) gdhA 8.085361e-06 0.0106215077 uvrA 1.730778e-01 0.7548230321 hdhA 1.064495e-08 0.0005746527 yecI 4.170175e-03 0.1651061600 galP 3.202789e-03 0.1471132276 malE 1.964593e-03 0.1187650888 Finally, to save the results, use 'write.csv': write.csv(results, file='myResults.csv', col.names=T, row.names=T) And you will have a csv file from your CyberT analysis written out to the directory you set earlier (for me, /Users/mkayala/develop/CyberT/).