RStudio character limit

I have been using my makebinary function without any issues, but now that I am working with a much larger tree, I am running into a problem with the number of characters in my newick tree.

My newick tree is 4434 characters long. Add to that the makebinary function call, and the whole command ends up being 4448 characters long. I found out that RStudio truncates commands to be 4096 characters. Fortunately, running it in the normal R console works without issues, but even that seems to have a 4096 character limit:

Note how there is now a new line (headed by the +) at the beginning of the second to last line.

R shortcuts

These might come up in a variety of contexts, but I think it's nice to know that R has a few built-in constants to prevent someone from tediously creating a vector such as c("a", "b", "c" ...). You can use LETTERS or letters, depending on whether you want the capital or lower-case letters. There are also abbreviated and full month names, month.abb and mouth.name, which are very handy for making plots. And finally, there is pi.

I recently used LETTERS and letters to create different plotting symbols for each individual data point so I could easily identify which point is what. I had more than 62 (26 capital letters, 26 lower-case letters, and 10 digits), so I had to use some symbols as well. I'm not sure if there's a more elegant way to do this, but if anyone knows of any, I'd love to hear about it.

for loops vs. sapply

Loops are very useful for doing the same (or similar things) multiple times. Unfortunately, in R, loops can be very clunky and slow.

For loops are perhaps more intuitive than sapply because the result you get is the same as if you ran the code within the for loop multiple times. What do I mean?

> x = numeric(10)

> y = numeric(10)

> z = for(i in 1:10) {

+ y[i] = i

+ x[i] = i*2

+ x[i]

+ }

> x

 [1]  2  4  6  8 10 12 14 16 18 20

> y

 [1]  1  2  3  4  5  6  7  8  9 10

> z

NULL

So the code within the for loop actually changes what is stored in x and y, but it does not return anything itself. Thus, z is NULL.

Let's use very similar code, except using sapply:

> x = numeric(10)

> y = numeric(10)

> z = sapply(1:10, function(i){

+    y[i] = i

+    x[i] = i*2

+    x[i]

+  })

> x

 [1] 0 0 0 0 0 0 0 0 0 0

> y

 [1] 0 0 0 0 0 0 0 0 0 0

> z

 [1]  2  4  6  8 10 12 14 16 18 20

Wait, why are x and y still full of 0s? This occurs because any assignments made within an sapply does not affect the global environment. So changing y[i] = i within sapply does not change the vector y itself. Thus it stays a vector of 0s as it was initialized. The trouble with sapply is that because of this, one iteration of the loop cannot depend on a different iteration of the loop--i.e., we cannot calculate x based off of what x was in a previous iteration. This is in direct contrast to a for loop, where because the changes happen in the global environment, we can use a previous iteration to determine the current iteration, like in this example:

> x = numeric(10)

> y = numeric(10)

> z = for(i in 2:10) {

+    y[i] = i

+    x[i] = x[i-1]+y[i-1]

+  }

> x

 [1]  0  0  2  5  9 14 20 27 35 44

> y

 [1]  0  2  3  4  5  6  7  8  9 10

> z

NULL

Fun with vectorizing

I found myself staying up way too late and having way too much fun helping a friend with her R homework. One of the things we did was to vectorize a function! Her assignment was to create plots using the Ricker model, with varying values for r:

So the first thing we did was to write a function that would calculate the population size over a certain time period for a single value of r:

Ricker = function(r, K = 100, n0 = 50, tf = 100) {
    n = numeric(length = tf)
    n[1] = n0
    for(t in 1:(tf-1)) {
          # We have tf - 1 because when t = 1, we are actually
          # calculating n[2], so if we want to end on n[tf], t 
          # should end on tf - 1
         n[t+1] = n[t]*exp(r*(1-n[t]/K))
    }
    n
}

Changing this to work with a longer vector of r values was fairly simple:

vectorizedRicker = function(r, K = 100, n0 = 50, tf = 100) {
     NumofRs = length(r)
      # First, we're changing our output to be
      # a matrix instead of a vector
     n = matrix(0, NumbofRs, tf)
      # We could vectorize over n0 as well, or any of the 
      # other parameters, if we set the number of rows for n
      # to be equal to the length of the longest vector
      # of parameter inputs
     n[,1] = n0
     for(t in 1:(tf-1)) {
           # Now, we just change the for loop to reflect that it
           # is working on entire columns of a matrix rather
           # than just a single element of a vector
          n[,t+1] = n[,t]*exp(r*(1 - n[,t]/K))
     }
     n
}

rfishbase

I was given a task to use the data in Fishbase to find percentages of different fish in different habitats. Carl Boettiger created an R package (rfishbase) to make the data very accessible in R, but unfortunately, some of the functions in the version on CRAN didn't work for me. Instead, I used the code he has on github. So if anyone else is having trouble getting updateCache() to work or wants to modify the code to extract a piece of information that isn't being extracted in the original fishbase() function, this is the place to go! Another nice thing is that you can also use getData() to download smaller chunks of data if you don't want to overwhelm the server.

Moving on from base

So I gave myself a project in my last post to remove hard-coded numbers in my code to plot a stacked histogram. I succeeded, but it may have been more work than it's worth. Perhaps it's time for me to move on to ggplot2 or lattice?

cats = levels(dat$Type) # same as before

xaxis = round(range(scores$x)*5)/5 # the range, rounded to every .2

breaks = seq(xaxis[1]-.2, xaxis[2]+.2, .2) # the breaks to use for histograms

histData = sapply(cats, function(x) {
     hist(scores$x[which(dat$Type == x)],
          breaks = breaks)$counts
})
allhistData = do.call(rbind, list(histData))
barplot(t(allhistData), space = 0, ylab = "number of sequences", xlab = "LD")
zero = which(breaks==0)-2 # find the 0 on my x-axis. The plot starts at -1

marks = seq(floor(xaxis[1]), ceiling(xaxis[2])) # the tick-marks I want to use, every integer

ticks = numeric() # this vector will hold where I want the tick labels to appear

while(zero >= -1) {
  ticks[length(ticks)+1] = zero
  zero = zero - 5
} # goes down by five, which corresponds to one LD unit

while(length(ticks) < length(marks)) {
  ticks[length(ticks)+1] = max(ticks) + 5
} # go up by five, until I have as many locations as I have labels

axis(side = 1, at = ticks[order(ticks)], labels = marks) #I need to order the tick marks in sequential order

legend("topleft", legend = cats, text.col = c("black", "gray"))

This produces the same plot from the previous post.

Linear Discriminant Analysis

This follows fairly naturally from the PCA I did on the data in my previous post. I have a dataset of several quantitative variables that can be grouped by a categorical variable. This time, I am going to maximize the separation between the two groups to see what traits are important in determining group inclusion.

While looking for information about linear discriminant analyses (LDA), I came across a very informative site by Dr. Avril Coghlan. I used a couple of the functions available from that website, including groupStandardise and calcWithinGroupsVariance, to obtain meaningful coefficients to determine what traits are more informative about group membership.

The lda() function is in the MASS package. Since I only have two groups, I only have one discriminant axis, as the number of discriminant axes is equal to the number of groups minus one. Thus I won't get a nice scatterplot the way I did with my PCA. Instead, if I try to plot the output to lda(), I get this:

I wanted a stacked histogram instead, so I had to do a little more work. 

cats = levels(dat$Type) # the categories I'm using

histData = sapply(cats, function(x) {

  hist(scores$x[which(dat$Type == x)],

       breaks = seq(-3, 2.2, .2))$counts

# bad me, I hard-coded the breaks. My next project can be to use the same breaks as the above plot

})

# this gives me a matrix of the counts for each interval for each type

allhistData = do.call(rbind, list(histData))

barplot(t(allhistData),

        space = 0,

        ylab = "number of sequences",

        xlab = "LD")

axis(side = 1,

     at = c(-1, 4, 9, 14, 19, 24),

     labels = c(-3, -2, -1, 0, 1, 2)) # more bad hard-coding

legend("topleft", legend = cats, text.col = c("black", "gray"))

And here is what I get:

It doesn't look great, but when did stacked histograms ever look good? Probably better if I have a lot more observations.

PCAs and Plotting

Principal Components Analysis, or PCA, is fairly straightforward using the princomp() function in R. But the data I have is divided by two factors: the type of response and the individual. I wanted to plot using different colors for the type of response and different plotting symbols for the individuals. Simple enough, but I also didn't want to need to modify the code if the number of types of responses or the number of individuals changed (the latter is probably more likely, but I still want my code to be as general as possible). Here is an example generated from random data:

The points themselves are fairly simple. I simply referred to the two variables I wanted in my col and pch arguments of plot().

col = as.numeric(dat[,x])+2

and

pch = dat$Individual+14

The reason for the +2 and +14 are to get the colors and plotting symbols I wanted. I could also assign specific colors I want by doing this:

colors = c("green", "blue")
ptype = c(15, 16, 17)

Then add these as arguments to plot().

col = colors[as.numeric(dat[,x])
pch = ptype[dat$Individual]

In this case, I would need to make sure that I have enough colors and plotting symbols that I don't run out.

The fun happens with the legend. For the text, I used this argument:

legend = c(levels(dat[,x]), unique(dat$Individual))

This remains flexible for any number of types or individuals. For the colors, I used a combination of seq() and rep() to get the numbers I wanted. If I didn't want my code to be general, I could simply use this:

text.col = c(3, 4, 1, 1, 1)

Instead, I used this:

text.col = c(seq(3, length(levels(dat[,x]))+2), rep(1, length(unique(dat$Individual))))

seq(3, length(levels(dat[,x]))+2) gives me a sequence of integers from three all the way to the number of types I have plus two (because I started with three instead of one). rep(1, length(unique(dat$Individual))) gives me 1 repeated for every unique individual.

Finally, the plotting symbols:

pch = c(rep(NA, length(levels(dat[,x]))), unique(dat$Individual)+14)

This is essentially the opposite of what I just did for the text color, except that I don't want any symbol next to the two types. Even though this looks like (and is) a lot more typing than simply hard-coding the appropriate numbers, this lets me use exactly the same code to make the figure even after I have doubled or tripled the amount of data I have.

RColorBrewer

Another great package that I was first introduced to in my Statistical Computing class, RColorBrewer.

This is a package that has built-in palettes that allows you to choose colors that have enough contrast for making plots.

The palettes come in three basic types: sequential, diverging, and qualitative. The sequential and diverging are great for plots where you want the colors to show an order. The difference between sequential and diverging seems to be a little subtle in terms of need: sequential shows more of a gradient, while diverging emphasizes both high and low extremes.

The qualitative palettes are best for categorical data with no ordering among categories. There are many sets, but they differ in the number of maximum colors, from 8 to 12. One interesting palette is the Paired palette, which consists of 6 hues, each with a light and dark color. I used this recently in a plot of different species, with males and females of varying lightness.

The same effect can be produced, perhaps to better effect, with different plotting symbols, but the Paired palette does a pretty good job.

likelihood reconstruction of ancestral states

In a previous post, I talked about how summarizing the state at each internal node over many make.simmap mappings did not correspond exactly with the ace reconstructions.  I contacted Liam Revell about this, and he informed me that this may be because ace does not compute the scaled marginal likelihoods but the conditional likelihoods of the subtrees descending from each node.  He suggested I try rerooting the tree at each internal node and using the ace reconstruction at the root to find the marginal likelihoods.  Here is my attempt:

nodes = (exampleTree$Nnode+2):(exampleTree$Nnode*2+1) # This gives me the number associated with each internal node

reRootAnc = t(sapply(nodes, function(x) {
  tr = reroot(exampleTree, node = x, position = 0) # rerooting the tree at each internal node
  reconst = ace(x = discrete, phy = tr, type = "discrete", model = "SYM") # estimating the maximum likelihood ancestral state estimate
  reconst$lik.anc[1,] # taking only the value for the root
}))

This is the result I got for an example tree:

Seems a little strange to me...I wonder what might be going on.

I knew this shouldn't be so complicated...

I was trying to analyze my incomplete dataset, and I needed to remove data for species where I have measured fewer than four individuals.  Since my dataset was small enough, it was easy to just remove them by hand, doing something like this:

> df
   Species meas
1        a    9
2        a    1
3        a    7
4        a    5
5        b    7
6        c    0
7        c    5
8        c    2
9        c    9
10       c    1
11       d    3
12       d    2

> dfremoved = df[-which(df$Species=="b"),]
> dfremoved = dfremoved[-which(dfremoved$Species=="d"),]
> dfremoved
   Species meas
1        a    9
2        a    1
3        a    7
4        a    5
6        c    0
7        c    5
8        c    2
9        c    9
10       c    1

But in order to do this systematically, I used a couple of steps.

> toofew = names(which(table(df$Species) < 4))
> toofew
[1] "b" "d"

First, I found the species names for species with fewer than four individuals.  With this, I can remove all rows where df$Species match any of these names.

> dfremoved = df[!(df$Species %in% toofew),]
> dfremoved
   Species meas
1        a    9
2        a    1
3        a    7
4        a    5
6        c    0
7        c    5
8        c    2
9        c    9
10       c    1

Avoiding Repetition

If there is anything I learned in STA 141, I learned the importance of the DRY principle:  don't repeat yourself.  Anything repetitive was heavily penalized in the grades, but it's also more prone to error and often takes much longer.  Even still, it's easy to be lazy and fall into the "copy/paste then change one word" routine, especially when just exploring a data set.  That's what I started out doing, but it turns out that doing it the 'right' way even easier!

For example, I want a plot of all morphological variables against size (standard length).  I want the points colored by species, and I want each point to be a unique value for the species.

plot(dat$standard.length, dat$head.length, col = as.factor(dat$Species), pch = as.character(dat$Number))

Now I can copy/paste this line and replace "head.length" with all of my other variables.  Simple enough.

But it turns out I have 24 variables.  So doing this will take much longer than this simple loop:

meas = names(dat)[9:length(dat)] # all morphological variables except standard length

sapply(meas, function(x) {
  png(file = paste(x, ".png", sep = ""))
  plot(dat$standard.length, dat[,x], col = as.factor(dat$Species), pch = as.character(dat$Number))
  dev.off()
})

I can do better by putting down axis labels and such, but now I have .png files of each of my morphological variables that I can browse through with my favorite image viewer.

Control Flow

I use if/else and for quite often, but rarely use while or repeat.  Even still, there are a couple of things with if/else that give me trouble if I haven't used it in a while.

If everything is on one line, all is good:

> x = 1

> y = if(x>0) 1 else 0

> y

[1] 1

But if I put the curly braces in the wrong place:

> if(x>0) {
+   x = -x
+   y = 1
+ }
> else {
Error: unexpected 'else' in "else"
>   x = x
>   y = 0
> }
Error: unexpected '}' in "}"

So I have to always make sure to do this:

> if(x>0) {

+   x = -x

+   y = 1

+ } else {

+   x = x

+   y = 0

+ }

> x

[1] -1

> y

[1] 1

No polytomies allowed?

I have recently been in a position to want an efficient way to obtain a binary Newick-format tree from a Newick-format tree that may or may not have polytomies.  Using a few functions from ape, this was fairly simple to obtain:


makebinary = function(newick) {
  tree = read.tree(text = newick)
  if(is.binary.tree(tree)) {
    return(newick)
  } else {
    return(write.tree(multi2di(tree)))
  }
}


This contains three very useful functions from ape:  read.tree, write.tree, is.binary.tree, and multi2di.  The first two are used to read/write Newick-format trees.  The functions read.nexus and write.nexus can be used for Nexus-format trees.  read.tree can be used for files as well.  is.binary.tree checks whether a tree has any polytomies, and multi2di converts a tree with multichotomies to a fully dichotomous tree with some branches of length 0.  There is also a function di2multi that will collapse any branches less than a tolerance level to a polytomy.

Animating SIMMAP trees

It took me quite a while to get this working, but I finally did:

The tree topology is from Rüber et al. (2004).  This is what I used to generate the image:

animateSimmap = function(phy, interval = .02, numb = length(phy), name = "animation.gif", ...) {
  if(class(phy) != "multiPhylo") stop("object must be multiPhylo")
  saveGIF(for(i in 1:numb) {
            dev.hold()
            plotSimmap(phy[[i]], ...)
            Sys.sleep(interval)
          }, movie.name = name)
}

So this function takes in a multiPhylo object that has stochastic character mappings, the time interval between trees (I'm not sure if this does anything because saveGIF might just have a set interval already), the number of mappings to animate, the name of the file to create, and any arguments you want to give to the plotSimmap function.  In terms of the mappings shown above, they're meaningless because the branches aren't proportional to anything meaningful, but it shows how this function would work.

plotSimmap is from the package phytools, and saveGIF is from the package animation.  To use saveGIF, you also need either ImageMagick (which is what I used), GraphicsMagick, or LyX.  

SIMMAP Trees

Thanks to Liam Revell, we can now produce simulated discrete character mappings within R.  The function make.simmap uses ace, from the package ape, to fit the model to use for the simulations.  Thus, the proportion at which a character state appears at a node among many iterations of simulations should be roughly equivalent to the likelihood of that state as estimated in ace.  I have heard this from many people, but I wanted to be able to summarize the actual simulated states.  In order to do this, I took a look at the components of a SIMMAP tree.  I'll illustrate this with a mock example:

We can take a look at what components this object has:

> names(exampleSimmap)
[1] "edge"        "edge.length" "tip.label"   "Nnode"       "maps"       
[6] "mapped.edge"

The first components are the same as any phylo object.  So maps and mapped.edge are what make a SIMMAP tree special.  Let's take a look (the middle elements removed to save space):

> exampleSimmap$maps
[[1]]
        1         0 
0.0958356 0.3197055 
[[2]]
        0         1 
0.2371619 0.1659281 
[[3]]
         1 
0.01380231 
...
[[17]]
        1 
0.1164989 
[[18]]
        1         0 
0.1161601 0.1769559 

This is precisely what we need!  Each element of exampleSimmap$maps represents a single branch, and the values represents the length of time that branch spends in each state, in this case 0 or 1.  That means we can simply take each branch's starting value (whatever is the name of the first element of that branch), and that is the value at the node where the branch starts.  Let's see if we can find this.

> exampleSimmap$edge
      [,1] [,2]
 [1,]   11   12
 [2,]   12   19
 [3,]   19    1
 [4,]   19    2
 [5,]   12   13
 [6,]   13   17
 [7,]   17    3
 [8,]   17    4
 [9,]   13   18
[10,]   18    5
[11,]   18    6
[12,]   11   14
[13,]   14   15
[14,]   15    7
[15,]   15   16
[16,]   16    8
[17,]   16    9
[18,]   14   10

The element named edge gives us the starting and ending node for each of the 18 edges in our tree.  That means we can use to figure out which node corresponds with which state.  Here I've written a function that takes in a SIMMAP tree and returns a named vector where the values are the node states and the names are the nodes.

mappedNode = function(phy) {
  # phy must be a SIMMAP tree
  nodes = phy$edge[,1] # this gives us the starting node for all edges
  map = sapply(phy$maps, function(x) attr(x, "names")[1]) # this gives us the starting value of each branch
  df = unique(data.frame(nodes = nodes, map = map)) # here we're removing the repeated values as interior nodes will have multiple branches
  mapping = df$map
  names(mapping) = df$nodes # naming the vector with node names
  mapping
}

Now we can run this function over all of the simulated mappings:

mappings = sapply(simmapTrees, function(x) mappedNode(x))
# change the 0 and 1 to numeric
nummaps = as.data.frame(sapply(1:length(simmapTrees), function(x) as.numeric(mappings[,x])))
# make sure the row names correspond to node names
rownames(nummaps) = rownames(mappings)
# get the number of simulated trees with  node state of 1
sums = sapply(rownames(nummaps), function(x) sum(nummaps[x,]))
# change that to a frequency
freq = sums/length(simmapTrees)
# plot
plot(exampleTree, label.offset = .05)
nodelabels(pie = freq, cex = .65, node = as.numeric(names(freq)))
tiplabels(pie = discTrait, cex = .65)
Looks fairly reasonable.  Now let's estimate the ancestral character using ace:

MLACE = ace(discTrait, exampleTree, type = "discrete", model = "SYM")

plot(exampleTree, label.offset = .05)
nodelabels(pie = 1 - MLACE$lik.anc, cex = .65)
tiplabels(pie = discTrait, cex = .65)

Pretty close, but not identical.  Maybe if we do more simulated mappings (I did 1000 here), they will start to look more similar?

Regular Expressions

My first experience with regular expressions came from Python for Dummies.  It looked particularly relevant to the specific task I was working on, scraping specific bits of information from fishbase.  When my advisor, Peter Wainwright, first approached me about this, I didn't know where to begin, so I went to two people with experience in these sorts of tasks:  Bob Thomson and Carl Boettiger.

Bob's suggestions were to download each individual fish's HTML file using a short bash script, then use something like Python or Perl to extract relevant bits.  With over 30,000 species, just downloading the HTML files took quite a long time.  But with no background (at the time) in Python or Perl, I turned to Carl, who suggested using R.  He quickly wrote a package, rfishbase that allows you to access information from the XML files on fishbase through R.  Although the XML files don't have all of the information available on the HTML files, they still have quite a lot.

My reason for this post, though, is because of a task my lab mate, Patrick, wished to accomplish using the data he accessed using rfishbase.  Looking at a character vector containing information of interest, he wanted to get all of the reference numbers within that vector.  An example of an element might be something like this:

"Occurs mainly over rocky and muddy bottoms.  Uncommon around coral reefs.  Usually rests on the bottom (Ref. 9710).  Juveniles may be found in shallow water, but adults are usually taken from depths of 70-330 m (Ref. 13442).  Reptant and natant decapods were the main food items throughout the year (Ref 59311).  Feeds on a wide variety of fishes and invertebrates."

Given this, he would want the numbers 9710, 13442, and 59311.  Even in this one example, you can see that they are not always consistent:  the first two have a period while the third doesn't.  And there are even things like this:

"Common species.  Free-living.  Assumed to feed on small invertebrates and fish (Ref.  4741, 34024).  Feed on small bottom animals (Ref. 35388)."

Notice the many spaces before the first reference and having two numbers.  Or this:

"Occurs in various inshore habitats (Ref. 9800). Feeds on benthic invertebrates and fish (Ref. 11889). Also Ref. 43081."

This one doesn't even have parentheses around the last one.  So the first thing I did was to find every instance of "Ref" followed by an optional period, any number of spaces, and a run of any number of numbers, commas, and spaces.

ref = regmatches(matches, gregexpr("Ref\\.? *[0-9 ,]*", matches))

where matches is the character vector with all of the information we're looking at.  This is modified from here.  Next, I removed all characters other than digits or commas and used strsplit to separate individual reference numbers.

refs = sapply(ref, function(x) unlist(strsplit(gsub("[^0-9,]", "", x), ",")))

You end up with a list the same length as the original character vector, and every element is a character vector of all of the reference numbers.  From here, you can go in and find all of the unique values to find all of the references you need.

Geiger Bug

Big thanks to Luke Harmon for supplying me with updated code that fixes a small bug in the geiger package!  For those of you who have had issues with fitDiscrete in the past, you can email Luke for a fix.

The issue that arises is when you use either the symmetrical model (model = "SYM") or the all rates different model (model = "ARD").  You will get an error message looking like this after waiting for some time for the likelihood optimization to finish:

Finding the maximum likelihood solution
[0        50      100]
[....................]
Error in getQ(exp(out$par), nRateCats, model) :
  You must supply the correct number of rate categories.

But you don't actually control what goes into getQ through the arguments you give to fitDiscrete.  So if you encounter this issue using either of the above models, be sure to email Luke to get the updated code.  A revamped version of geiger is on the way, so hopefully this issue won't be around for much longer!

The one thing Luke cautioned me about using these two models is that they can both quickly become parameter-rich.  The number of rate parameters for the symmetric model is n * (n-1) / 2, while for the all rates different model, the number of rate parameters is n * (n-1) where n is the number of discrete states.  For example, if I have five discrete states, I would have ten parameters for the symmetric model and twenty parameters for the all rates different model.  I need a lot of data to be estimating so many parameters!

I did try out his updated code, and it works perfectly fine.  It still does take a while (I'm sure the speed depends on the size of your data set, the shape of your phylogeny, your computer's specs...), but I hear that there is already a faster version of the code if you ask Luke for it.  It would definitely be interesting to test out how much of a difference there is and what they changed to make it faster.

If all you are trying to estimate is a single rate model (model = "ER", which stands for equal rates), then there is no need to use an updated version of the function.  The old version will work just fine.  The problem with the original was that the calculation of the rate categories occurred twice:  once to get the number of rate categories, and another within getQ to test whether the number of rate categories was correct.  So for example, with five discrete states nRateCats in the above will equal ten for the symmetric model.  But you will get an error to supply the correct number of rate categories because within getQ, it compares the number of parameters to nRateCats * (nRateCats-1) / 2, which would compare ten to forty-five.

Creating Trees

I told Luke Mahler about how I'm interested in practicing analyses in R with my particular group of interest, even though I haven't collected data on more than a few species and I don't have a tree.  He told me creating a tree from previously published studies is easy:  just import the tree into R in newick format! In the geiger package in R, there is a function read.tree that does just that.

The tricky part about doing this is creating the newick-format tree itself.  It involves a lot of parentheses, colons, and parentheses:

exampleTree = read.tree(text = "(((A:1, B:1):1, C:2):1, D:3);")

will give you

Since I was typing out the newick-format tree directly within read.tree, I needed to make sure it went into the text argument instead of the file argument.  As for the newick format itself, the tip names are followed by the immediately subtending branch length:  hence A:1 or D:3.  The nodes also need to be provided with branch lengths, which is the :1 following each clade designated in parentheses. Finally, you can't forget the semi-colon at the end.  It is very easy to get lost in a sea of parentheses, especially for large trees, so inevitably I needed to go back to make changes.  Testing small clades at a time makes this a little easier.

In my particular tree, I had polytomies, which I designated with 0-length internal branches.  So if in the above tree, I actually don't have any information about the interrelationships of A, B, and C, I can use this:

polytomyTree = read.tree(text = "(((A:1, B:1):0, C:1):1, D:2);")

to get

RStudio

I was first introduced to RStudio when I took STA141 with Duncan Temple Lang.  At the time, I was using a PC, and it made things much simpler.  A few of the features that I love about it:

1) Color scheme for scripts.  I can see at a glance what's a comment, and the parentheses pop out at you.  This was something I envied about Macs when I used to use my PC exclusively.  I do wish you could personalize it a little more, but the schemes they have are great.  I personally like Cobalt.

2) Balancing parentheses and quotations.  Yes, this can get obnoxious if you are tweaking existing code rather than starting from scratch, since if you try to type an end quote, it will be interpreted as another beginning quote.  But I find being able to see what's matching your current end parenthesis incredibly helpful.  Besides, if you don't like it, you can easily turn it off in the preferences.

3) Everything in one window.  This wasn't a step up from my PC version since the normal R console keeps everything in one window, but it's so nice to have everything (script, console, history, help pages, graphics) in one window.