Published on By Tim Sackton , tagged as R

R is a functional programming language, which means that most of what one does is apply functions to objects.

We will begin with a brief introduction to R objects and how functions work, and then focus on getting data into R, manipulating that data in R, and generating basic summary statistics.

## Getting Help

R has very good built-in documentation that describes what functions do.

To get help about a particular function, use ? followed by the function name, like so:

```?read.table
```

## Vectors

Let’s start by creating one of the simplest R objects, a vector of numbers.

```v1<-c(1,2,3,4,5)
```

v1 is an object which we created by using the `<-` operator (less followed by a dash).

v1 now contains the output of the function `c(1,2,3,4,5)` (c, for combined, combines elements into a vector

Let’s display the contents of v1:

```print(v1)

##  1 2 3 4 5
```

We can also just type the name of an object (in this case `v1`) to display it.

Let’s make a few more objects:

```x<-10
y<-11
letters<-c("a", "b", "c", "d")
```

We might want to get a list of all the objects we’ve created at this point. In R, the function ls() returns a character vector with the names of all the objects in a specified environment (by default, the set of user-defined objects and functions).

```ls()

##  "letters" "v1"      "x"       "y"
```

We can manipulate objects like so:

```x+5

##  15

print(x)

##  10

x*2

##  20

print(x)

##  10
```

Note that the value of `x` is not modified here. We did not save the output to a new object, so it is printed to the screen.

If we want to update the value of `x`, we need to use our assignment operator:

```x<-x+5
print(x)

##  15
```

R handles vector math for us automatically:

```v1*x

##  15 30 45 60 75

v2<-v1*x
print(v2)

##  15 30 45 60 75
```

## Object Types

All objects have a type. Object types are a complex topic and we are only going to scratch the surface today. To slightly simplify, all data objects in R are either atomic vectors (contain only a single type of data), or data structures that combine atomic vectors in various ways. We’ll walk through a few examples. First, let’s consider a couple of the vectors that we’ve already made.

```class(x)

##  "numeric"

class(letters)

##  "character"
```

Numeric and character data types are two of the most common we’ll encounter, and are just what they sound like. Another useful type is logical data, as in TRUE/FALSE. We can create a logical vector directly like so:

```logic1<-c(TRUE, TRUE, FALSE, FALSE)
print(logic1)

##   TRUE  TRUE FALSE FALSE
```

Or we can use logical tests.

```logic2<-v1>2
print(logic2)

##  FALSE FALSE  TRUE  TRUE  TRUE
```

Note that the logical test here (>2) is applied independently to each element in the vector. We’ll come back to this when we talk about data subsets.

A final thing to note about logical vectors: they can be converted to numeric with TRUE=1 and FALSE=0. So a simple way to find how many elements in a logical vector are true is `sum()`.

```sum(logic2)

##  3
```

## Data Frames

So far we’ve been talking about atomic vectors, which only contain a single data type (every element is logical, or character, or numeric). However, data sets will usually have multiple different data types: numeric for continunous data, character for categorical data and sample labels. Depending on how underlying types are combined, we can have four different “higher-level” data types in R:

Dimensions Homogeneous Heterogeneous
1-D atomic vector list
2-D matrix data frame

We’ll focus on data frames, as those are the data type most commonly found in real biological analysis. A data frame is a collection of vectors, which can be (but don’t have to be) different types, but all have to have the same length. Let’s make a couple of toy data frames.

One way to do this is with the data.frame() function.

```df1<-data.frame(label=c("rep1", "rep2", "rep3", "rep4"), data=c(23, 34, 15, 19))
class(df1)

##  "data.frame"

print(df1)

##   label data
## 1  rep1   23
## 2  rep2   34
## 3  rep3   15
## 4  rep4   19
```

str() gives lots of information about the data type of the consituent parts of a data frame

```str(df1)

## 'data.frame':    4 obs. of  2 variables:
##  \$ label: Factor w/ 4 levels "rep1","rep2",..: 1 2 3 4
##  \$ data : num  23 34 15 19
```

Note that label is a factor – this is a special kind of character vector. By default, when creating a data frame or reading data from a file, R will convert strings (character vectors) to factors. Factors can be useful, but also introduce some pitfalls, which we can come back to if we have time. A good source of information on factors is here.

A more common way to create a data frame is by reading from a file. There are a few functions to do this in R: read.table() and read.csv() are some of the most common.

```bus<-read.table(file="http://software.rc.fas.harvard.edu/ngsdata/workshops/2015_March/mbta_bus.tsv", header=TRUE, sep="\t", stringsAsFactors=T)

#note that R can read data directly from the web -- no need to download the file first

##   route      type cost.per.pax ridership pax.per.trip
## 1     1       Key         0.63     13306          416
## 2     4  Commuter         4.38       459           21
## 3     5 Community         2.99       156            5
## 4     7     Local         1.52      3893          177
## 5     8     Local         2.02      3844           85
## 6     9     Local         1.21      5980          187
```

`head()` displays the first few rows of a data frame and is great for examining large objects. Let’s also look at `str()`

```str(bus)

## 'data.frame':    163 obs. of  5 variables:
##  \$ route       : Factor w/ 163 levels "1","10","100",..: 1 74 104 128 142 151 2 10 25 26 ...
##  \$ type        : Factor w/ 5 levels "Community","Commuter",..: 4 2 1 5 5 5 5 5 5 4 ...
##  \$ cost.per.pax: num  0.63 4.38 2.99 1.52 2.02 1.21 1.75 1.84 2.66 1.11 ...
##  \$ ridership   : int  13306 459 156 3893 3844 5980 3184 3312 1285 6227 ...
##  \$ pax.per.trip: int  416 21 5 177 85 187 398 237 22 111 ...
```

summary() is a good way to see some basic information about a data frame

```summary(bus)

##      route            type      cost.per.pax     ridership
##  1      :  1   Community:  1   Min.   :0.570   Min.   :   30.0
##  10     :  1   Commuter : 13   1st Qu.:1.350   1st Qu.:  755.5
##  100    :  1   Express  : 10   Median :2.150   Median : 1481.0
##  101    :  1   Key      : 15   Mean   :2.433   Mean   : 2427.1
##  104    :  1   Local    :124   3rd Qu.:3.055   3rd Qu.: 3101.0
##  105    :  1                   Max.   :9.760   Max.   :15018.0
##  (Other):157
##   pax.per.trip
##  Min.   :   0.00
##  1st Qu.:   9.00
##  Median :  21.00
##  Mean   :  78.55
##  3rd Qu.:  48.00
##  Max.   :2550.00
##
```

## Subsetting Data

Now let’s look at ways to get subsets of data. R uses `[]` to index rows and columns, like so:

```bus[1,]

##   route type cost.per.pax ridership pax.per.trip
## 1     1  Key         0.63     13306          416

##  1 4 5 7 8 9
## 163 Levels: 1 10 100 101 104 105 106 108 109 11 110 111 112 114 116 ... CT3

bus[1,1]

##  1
## 163 Levels: 1 10 100 101 104 105 106 108 109 11 110 111 112 114 116 ... CT3

bus[15,2]

##  Local
## Levels: Community Commuter Express Key Local

bus[1:5,4]

##  13306   459   156  3893  3844

bus[c(1,3,5),c(2,3)]

##        type cost.per.pax
## 1       Key         0.63
## 3 Community         2.99
## 5     Local         2.02
```

Of course, usually we want to select a set of rows from a data frame matching some criteria. We can use the subset() function for this. Let’s start by just getting the data for key bus routes.

```key.bus<-subset(bus, type=="Key")
```

Sometimes we might want something more complicated. Let’s say we want to get data for all “local” bus routes with at least 1000 daily riders.

```local.bus<-subset(bus, type=="Local" & ridership > 1000)
```

We can also use [] (see `?Extract`) to subset data in more complicated ways. The key insight here is that we can use logical vectors as indexes. So for example:

```high.ridership<-bus\$ridership > 5000
class(high.ridership)

##  "logical"

print(high.ridership)

##     TRUE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE  TRUE
##   FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE  TRUE FALSE
##   FALSE FALSE  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE FALSE
##   FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
##   FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE  TRUE
##   FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE
##   FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##   FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##    TRUE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

bus[high.ridership,]

##    route  type cost.per.pax ridership pax.per.trip
## 1      1   Key         0.63     13306          416
## 6      9 Local         1.21      5980          187
## 10    15   Key         1.11      6227          111
## 11    16 Local         0.99      5100         2550
## 16    22   Key         1.05      8151          170
## 17    23   Key         1.00     11687          899
## 21    28   Key         0.80     13550          387
## 25    32   Key         1.02     10130          103
## 27    34 Local         2.43      5959           60
## 32    39   Key         0.72     15018          406
## 44    57   Key         0.86     12081          263
## 50    66   Key         0.79     13630          147
## 54    70 Local         1.96      7021          195
## 55    71   Key         1.35      5488          157
## 57    73   Key         1.75      5834           46
## 61    77   Key         1.95      7253          165
## 68    86 Local         1.07      5844           76
## 89   111   Key         1.06     11695           38
## 92   116   Key         0.88      5077           26
```

We don’t need to store the logical vector in an object, so we can get the same result with this:

```bus[bus\$ridership > 5000,]

##    route  type cost.per.pax ridership pax.per.trip
## 1      1   Key         0.63     13306          416
## 6      9 Local         1.21      5980          187
## 10    15   Key         1.11      6227          111
## 11    16 Local         0.99      5100         2550
## 16    22   Key         1.05      8151          170
## 17    23   Key         1.00     11687          899
## 21    28   Key         0.80     13550          387
## 25    32   Key         1.02     10130          103
## 27    34 Local         2.43      5959           60
## 32    39   Key         0.72     15018          406
## 44    57   Key         0.86     12081          263
## 50    66   Key         0.79     13630          147
## 54    70 Local         1.96      7021          195
## 55    71   Key         1.35      5488          157
## 57    73   Key         1.75      5834           46
## 61    77   Key         1.95      7253          165
## 68    86 Local         1.07      5844           76
## 89   111   Key         1.06     11695           38
## 92   116   Key         0.88      5077           26
```

We might want to add some new variables to our dataset. For example, perhaps we want to know the cost per day to operate each bus line. We can calculate this from the ridership and cost.per.pax numbers.

```bus\$cost.per.day<-bus\$cost.per.pax * bus\$ridership

##   route      type cost.per.pax ridership pax.per.trip cost.per.day
## 1     1       Key         0.63     13306          416      8382.78
## 2     4  Commuter         4.38       459           21      2010.42
## 3     5 Community         2.99       156            5       466.44
## 4     7     Local         1.52      3893          177      5917.36
## 5     8     Local         2.02      3844           85      7764.88
## 6     9     Local         1.21      5980          187      7235.80
```

Let’s also compute the revenue per day for each bus line. This is actually a bit complicated because not all buses charge the same rate. For local/key/community lines, it is \$1.60, but for express/commuter lines it is either \$3.65 or \$5.25. So we do this in a few steps. Let’s first make a fare column that lists the fare for each line.

```bus\$fare<-ifelse(bus\$type != "Commuter" & bus\$type != "Express", 1.60, 3.65)
bus\$fare[bus\$route=="352" | bus\$route=="354" | bus\$route=="505"] = 5.25
```

Now we can calculate revenue per line as fare * ridership.

```bus\$revenue.per.day = bus\$fare * bus\$ridership
```

Now let’s explore this dataset a bit. What is the mean ridership on MBTA bus lines?

```mean(bus\$ridership)

##  2427.055
```

What if we wanted to get the means for all the (numeric) columns? We could type each by hand, but R has a shortcut – the apply function.

```apply(bus[,c(3,4,5,6)], MARGIN=2, mean)

## cost.per.pax    ridership pax.per.trip cost.per.day
##     2.433067  2427.055215    78.552147  3755.512883
```

Apply can be used to apply a function (in this case mean) to all the columns (MARGIN=2) or to all the rows (MARGIN=1)

## Writing Data

We’ve added several variables now, and we might want to write our updated dataset to a file. We do this with the write.table() function, which writes out a data.frame.

```write.table(bus, file="mbta_bus_updated.tsv", quote=F, sep="\t", row.names=F, col.names=T)
```

But where did R put the file on our computer? R does all file operations without a full path in a working directory. RStudio has a preference to set the default working directory, which is typically something like `/Users/Tim/R`.

To see the current working directory, use:

```getwd()

##  "/Users/williams03/a/workshops/2015.March/Intro_to_R"
```

We can change the working directory with `setwd()`.

## Plotting

Now, let’s do some visualizations. We’ll start with something simple – making a histogram.

How about the distribution? We’ll use the hist() function to plot a histogram

```hist(bus\$ridership, col="gray50", breaks=15, xlab="Ridership", main="")
```

Another very useful way to summarize data is with a boxplot. Here, we’ll look at ridership for each bus type. However, we’ll ignore the one bus line with type==Community.

```boxplot(ridership ~ type, data=bus[bus\$type!="Community",])
```

Note that Community still shows up on the axis, even though we are not plotting the data for that bus line. This is because ‘type’ is stored as a factor, and all levels of a factor (even those with no actual data) are plotted. One simple way to remove factor levels that don’t appear in the dataset is with `droplevels()`.

```boxplot(ridership ~ type, data=droplevels(bus[bus\$type!="Community",]))
```

Finally, let’s look at an X-Y plot, in this case comparing ridership to cost per passenger.

```plot(x=bus\$ridership, y=bus\$cost.per.pax, pch=16)
```

Looks like a non-linear relationship, so let’s try a log-log plot

```plot(x=bus\$ridership, y=bus\$cost.per.pax, pch=16, log="xy")
```

Alternatively:

```plot(log10(bus\$cost.per.pax) ~ log10(bus\$ridership), pch=16)
```

## Statistics

Now let’s demonstrate a few statistical tests.

Our boxplot suggests that ‘key’ bus routes have much higher ridership than other bus routes. Let’s test that using a T-test.

```#note that we can use a logical test as a grouping variable in a formula
t.test(ridership ~ type=="Key", data=bus)

##
##  Welch Two Sample t-test
##
## data:  ridership by type == "Key"
## t = -8.3121, df = 14.435, p-value = 7.043e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -9905.928 -5851.498
## sample estimates:
## mean in group FALSE  mean in group TRUE
##            1702.020            9580.733
```

We might want to do this non-parameterically using a Wilcoxon Rank Sum Test.

```wilcox.test(ridership ~ type=="Key", data=bus)

##
##  Wilcoxon rank sum test with continuity correction
##
## data:  ridership by type == "Key"
## W = 22, p-value = 4.281e-10
## alternative hypothesis: true location shift is not equal to 0
```

Our plot suggested that log(ridership) and log(cost per passenger) are strongly correlated. Let’s test that with cor.test()

```cor.test(log10(bus\$ridership), log10(bus\$cost.per.pax))

##
##  Pearson's product-moment correlation
##
## data:  log10(bus\$ridership) and log10(bus\$cost.per.pax)
## t = -13.3327, df = 161, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7901266 -0.6421777
## sample estimates:
##        cor
## -0.7243894
```

Now let’s work through a final example that involves both data manipulation and some more statistical tests. Let’s say we want to find out how many bus lines are profitable (make more in fare revenue than they cost to operate)? What would you do?

```bus\$cost.per.day < bus\$revenue.per.day

##     TRUE FALSE FALSE  TRUE FALSE  TRUE FALSE FALSE FALSE  TRUE  TRUE
##    TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
##    TRUE  TRUE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE  TRUE  TRUE
##   FALSE  TRUE FALSE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE
##   FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE  TRUE FALSE  TRUE
##   FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE
##   FALSE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE FALSE  TRUE FALSE FALSE
##   FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
##    TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE FALSE FALSE FALSE
##  FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE
##  FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  FALSE FALSE  TRUE  TRUE FALSE  TRUE FALSE  TRUE FALSE FALSE FALSE
##   TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
##  FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE
##  FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE

#this is a logical vector; remember sum() will give us the number of elements that are TRUE
sum(bus\$cost.per.day < bus\$revenue.per.day)

##  63
```

Now we might want to know if different types of buses are more or less likely to be profitable. Let’s start by making a new column in our data frame that is TRUE if a bus route is profitable and FALSE otherwise. How would you do that?

```bus\$profitable = bus\$cost.per.day < bus\$revenue.per.day
```

Now we can use the table() function to look for associations between categorical variables.

```table(bus\$profitable, bus\$type)

##
##         Community Commuter Express Key Local
##   FALSE         1        8       2   2    87
##   TRUE          0        5       8  13    37

#to test non-independence using an RxC chisquare test, use chisq.test
chisq.test(table(bus\$profitable, bus\$type))

## Warning in chisq.test(table(bus\$profitable, bus\$type)): Chi-squared
## approximation may be incorrect

##
##  Pearson's Chi-squared test
##
## data:  table(bus\$profitable, bus\$type)
## X-squared = 26.4862, df = 4, p-value = 2.525e-05
```

Notice we got a warning message here. In this case, this means that our data violates the assumptions of the standard Chi-squared approximation (likely because there are several cells with low counts). In general, it is always wise to understand the cause of any warnings you get.

Finally, we’ll end with a note about missing data. R stores missing data with the special value `NA`. While our toy dataset here has no missing values, this will usually not be the case with real data. You can read more about missing data and data types here

## Preparation for next time

For next Tuesday, we will be using packages from Bioconductor. Packages extend R in various ways – Bioconductor, in particular, contains a very useful set of packages for handling various kinds of genomic and high-throughput biological data.

The easiest way to install packages in RStudio is using the pull-down menu, but you can also use the function `install.packages()`. For example, `install.packages(“maps”)` will install a package for using R to draw maps.

Packages, however, are not automatically loaded. To load a package, use `library(packagename`). For example,

```library(maps) #loads R package to draw geographic maps
```

Installing packages from Bioconductor is a little more complicated and requires running a helper script from the Bioconductor web site, but this is actually pretty easy:

```source("http://bioconductor.org/biocLite.R")
```

Now we can use the `biocLite()` function to install Bioconductor packages.

Important: prior to next Tuesday, please install DESeq2 and edgeR packages from Bioconductor, with these commands:

```biocLite("DESeq2")
biocLite("edgeR")
```

This can take several minutes to run as if you have never installed any Bioconductor packages you will have to install a number of prerequisites.