So you want to use your computer for science…

Posted on November 5, 2017 by Jeff

It’s been a while since I was a new graduate student, and I’ve forgotten how little I knew about computers back then.  I was reminded recently while teaching a couple of lab members how to use ffmpeg, an excellent command line tool for building animations from images (as described in this post).  We got there, but I realized that we needed a basic computing tutorial before moving on to anything more advanced.  If you’re trying to use your laptop for science, but you’re not too sure about this whole “command line thing”, this post’s for you.  Be warned that this tutorial is intended as only the most cursory crash course to get you moving up the initial learning curve.  A comprehensive list of commands for Bash on OSX can be found here.  At the end of this tutorial you will have a basic grasp of how to:

  • Navigate the file system
  • Create and edit text files
  • Search text files
  • Create a shell script
  • Modify the PATH variable using startup files (.bash_profile)

Okay, so I want to get a computer to do some science.  What should I get?

Whatever you want.  It really doesn’t matter.  Most grad students seem to get Macs these days, which I don’t love (they’re costly and epitomize form over function), but have the slight advantage of a Unix-like environment hiding behind all that OSX gibberish.  I use a Windows machine (which I also don’t love, as it epitomizes dysfunction over function) with Cygwin, which gives me access to all the Linux tools that I need to carry out day-to-day operations.  Windows 10 users can also make use of the Bash shell add-on for Windows, but I haven’t found any advantage to this over Cygwin.  The point is that you need either a) A Mac, b) A Windows machine running Cygwin or the add-on, or c) A Linux machine.  The command prompt and output given below are what you will see in OSX (faked since I’m not actually using a Mac), but are similar to what you will get with Cygwin or one of the Linux distros.  The commands should work the same across all of these options.

Getting familiar with “Bash”

Bash stands for Bourne-again shell, which you can read all about in many other places on the web.  Bash is a very powerful tool for manipulating your computer’s file system, executing programs, and even creating programs.  It is a cornerstone of scientific computing and you should have at least some passing familiarity with it.  To open the Bash terminal in OSX, go to Applications/Utilities/Terminal in Finder.  A mysterious black (or white) window will open, with a white (black) cursor waiting for YOU.  Type “pwd” for print working directory and hit “Enter”.

jeffscomputer:~ jeff$ pwd
/Users/jeff

Bash will respond with your location, which should be something along the lines of /user/home.  Next type “ls” to list the contents of the directory.

jeffscomputer:~ jeff$ ls
bin
Desktop
Documents
Downloads

We can use the command “cd” to change directories.  For example, if you want to move to your Desktop directory type “cd Desktop”.

jeffscomputer:~ jeff$ cd Desktop
jeffscomputer:Desktop jeff$ pwd
/Users/jeff/Desktop

Because the directory “Desktop” was just one level down, in this case the relative path “Desktop” is equivalent to typing the full path “/Users/jeff/Desktop”.  Here’s a useful tip.  From any location on your system “cd ~” will get you back to your home directory.

jeffscomputer:Desktop jeff$ cd ~
jeffscomputer:~ jeff$ pwd
/Users/jeff

Now let’s create a directory to do some work.  The command for this is “mkdir temp”, for “make directory with name temp.”

jeffscomputer:~ jeff$ mkdir temp
jeffscomputer:~ jeff$ cd temp

Now move into that directory.  You already know how 🙂

Creating and editing text files

You will frequently need to create and edit basic text files without all the fancy formatting of a word processing document.  The most user friendly way to do this is with the program nano, which is likely already present if you are using OSX or Cygwin.  Type “nano temp.txt” and nano will open a blank text file with name temp.txt.

jeffscomputer:temp jeff$ nano temp.txt

Type a couple lines of text and, when you’re ready to exit and save the file, hit ctrl-x.  Nano will prompt you about saving the output, hit yes.  List the contents of the directory and notice that the file temp.txt now exists.  Type “nano temp.txt” again and rather than create a new file, nano will open temp.txt for editing.

Having gone through the trouble of creating that file, let’s go ahead and delete it using the remove or “rm” command.

jeffscomputer:temp jeff$ rm temp.txt

To do some fancier things with files lets download one that has a little more information in it.  There are two programs to use to fetch files online, wget and curl.  I find wget much easier to use than curl.  If you’re using Cygwin or Linux you likely already have it, but for OSX you first need to install a package manager, which is a whole can of worms that I’m not going to tackle in this post.  So let’s use curl to download a text file, in this case “The Rime of the Ancient Mariner” by Samuel Coleridge.  The basic download command for curl is:

jeffscomputer:temp jeff$ curl -O https://www.polarmicrobes.org/extras/ancient_mariner.txt

This should create the file “ancient_mariner.txt” in your working directory (e.g., “/Users/jeff/temp”).

Viewing file content

The reason we downloaded this more complex text file (it’s a pretty long poem) is to simulate a longer data file than our “temp.txt” file.  Very often in scientific computing you have text files with hundreds, thousands, or even millions of lines.  Just opening such a file can be onerous, let alone finding some specific piece of information.  Fortunately there are tools that can help.  Type “head ancient_mariner.txt”, this returns the top 10 lines of the file.

jeffscomputer:temp jeff$ head ancient_mariner.txt
PART THE FIRST
It is an ancient mariner,
And he stoppeth one of three.
"By thy long grey beard and glittering eye,
Now wherefore stopp'st thou me?
"The Bridegroom's doors are opened wide,
And I am next of kin;
The guests are met, the feast is set:
May'st hear the merry din."
He holds him with his skinny hand,

Want to guess what “tail” does?  For either command you can use a “flag” to modify behavior, such as returning more lines.  Flags are always preceded by “-” or “–“, and generally come before the positional arguments of the command, in this case the file.  This is general syntax for Unix commands, and does not apply only to “head” and “tail”.

jeffscomputer:temp jeff$ tail -15 ancient_mariner.txt
Both man and bird and beast.
He prayeth best, who loveth best
All things both great and small;
For the dear God who loveth us,
He made and loveth all.
The Mariner, whose eye is bright,
Whose beard with age is hoar,
Is gone: and now the Wedding-Guest
Turned from the bridegroom's door.
He went like one that hath been stunned,
And is of sense forlorn:
A sadder and a wiser man,
He rose the morrow morn.

In the examples so far the the command output has printed to the screen, but what if we want to capture it in a file?  For example, what if we have a huge data file (millions of rows) and we want just the top few lines to test some code or share with a collaborator?  This is easily done by redirecting the output using “>”.

jeffscomputer:temp jeff$ tail -15 ancient_mariner.txt > end_of_ancient_mariner.txt

Searching file content

To find specific information in a file use the command “grep”.  Without launching into a full-on explanation of regular expressions, grep (very quickly) finds lines that match some given pattern.  You can either count or view the lines that match.  To find the lines with the word “albatross” in ancient_mariner.txt:

jeffscomputer:temp jeff$ grep 'Albatross' ancient_mariner.txt
At length did cross and Albatross,
I shot the Albatross."
Instead of the cross, the Albatross
The Albatross fell off, and sank
The harmless Albatross."
The Albatross's blood.

Seems there’s a typo in this version of the poem (and|an)!  At any rate, use the -c flag to count the lines.  Protip: use the “up” key to bring back the previous line, which you can then modify.

jeffscomputer:temp jeff$ grep -c 'Albatross' ancient_mariner.txt
6

Use the -v flag to select against lines with “Albatross”, which you can combine with the -c or other flags:

jeffscomputer:temp jeff$ grep -cv 'Albatross' ancient_mariner.txt
634

Build a basic program

So far we’ve been executing commands manually, from the command line.  Suppose we have a set of commands that we want to execute a lot, or that need a method to document our workflow.  To do this we create a shell script.  Fire up nano for a file named “temp.sh” and type:

#!/bin/bash

echo "hello world!"
# hey, this line doesn't do anything!
f=ancient_mariner.txt
echo $f
grep -cv 'Albatross' $f

Line by line here’s what’s going on:

  • The first line is called the shebang and it tells your system what interpreter to use to run the script (in this case Bash).  /bin/bash is an actual location on your computer where the Bash program resides
  • The next line is just a bit of formality; by strictly adhered-to convention your first program should always be a little script that says “hello world!”.  It does however, illustrate the “echo” command, which prints out information.
  • The next line starts with “#” which denotes a comment.  Line that start with “#” are not read by Bash.  You can use that character to make notes, or to toggle commands on and off.
  • In the next line we assign a variable (f) a value (ancient_mariner.txt).  We can now call that variable using “$”.  The next two lines are examples of this.

To execute the script we simply type the file name into bash.  Before we do that however, we need to set the file permissions, as files are not by fault executable.  To do that we use the “chmod” command with the “a+x” options (note that this is not a flag).

jeffscomputer:temp jeff$ chmod a+x temp.sh

You can run it now, but there is one final trick.  Bash doesn’t know to look in your working directory for the script, you have to specify that that’s where it is.  The location of the current working directory is always “./”, so the command looks like this:

jeffscomputer:temp jeff$ ./temp.sh
hello world!
ancient_mariner.txt
634

Setting up your environment

Okay, we’re going to ramp things up a bit for the grand finale and modify the Bash startup files to better use your new-found skills.  There are several possible startup files, and the whole startup file situation gets a bit confusing.   We will modify the .bash_profile file, which will handle the majority of user cases, but you should take the time to familiarize yourself with the different files.  The .bash_profile file is a hidden file (denoted by the .), you can see hidden files by using “ls” with the “-a” flag.

jeffscomputer:temp jeff$ cd ~
jeffscomputer:~ jeff$ ls -a
bin
Desktop
Documents
Downloads
.bash_history
.bashrc
.profile
.ssh

Don’t worry if you don’t see .bash_profile, we will create it shortly.  First, to understand why we need to modify the startup file in the first place, from your home directory try executing the shell script that we just created.

jeffscomputer:~ jeff$ temp.sh
-bash: temp.sh: command not found

The command would execute if you typed temp/temp.sh, but remembering the location of every script and program so that you can specify the complete path to it would be silly.  Instead, Bash stores this information in a variable named PATH.  To see what’s in PATH use “echo”:

jeffscomputer:~ jeff$ echo $PATH
/usr/local/bin:/usr/bin:/usr/sbin

If you created a script in /usr/local/bin, Bash would know to look there and would find the script.  It’s a good idea however, to keep user-generated scripts in the home directory to avoid cluttering up locations used by the operating system.  What we need to do is automatically update PATH with customized locations whenever we start a new Bash session.  We accomplish this by modifying PATH on startup using the startup file.  To do this use nano, but you will need to use nano as a super user or “sudo”.

jeffscomputer:~ jeff$ sudo nano .bash_profile

As you already know, if you don’t have .bash_profile this will create it for you.  Now, at the end of the .bash_profile file (or on the first line, if its empty), type (replacing “jeff” with your home directory):

export PATH=/Users/jeff/temp:$PATH

The command structure might look opaque at first but its really not.  This line is saying “export the variable PATH as this text, followed by the original PATH variable”.

Close nano.  Recall that .bash_profile is read by Bash at the start of the Bash session.  Your newly defined PATH variable will be read if you start a new session, but you can also force Bash to read it with the “source” command.

jeffscomputer:~ jeff$ source .bash_profile

Now try to execute your bash script.

jeffscomputer:~ jeff$ temp.sh
hello world!
ancient_mariner.txt
634

Last, let’s clean things up a little.  Previously we used “rm” to remove a file.  The same command with the -r flag will remove a directory.

jeffscomputer:~ jeff$ rm -r temp

Voila!  To keep your .bash_profile file looking pretty, don’t forget to remove the line adding temp to PATH (though it does no harm).  Or you can comment that line out and leave it as an example of the correct syntax for when you next add a new location to PATH.

Posted in Computer tutorials | Leave a comment

Analyzing flow cytometry data with R

Posted on August 11, 2017 by Jeff

We recently got our CyFlow Space flow cytometer in the lab and have been working out the kinks.  From a flow cytometry perspective the California coastal environment is pretty different from the western Antarctic Peninsula where I’ve done most of my flow cytometry work.  Getting my eyes calibrated to a new flow cytometer and a the coastal California environment has been an experience.  Helping me on this task is Tia Rabsatt, a SURF REU student from the US Virgin Islands.  Tia will be heading home in a couple of weeks which presents a challenge; once she leaves she won’t have access to the proprietary software that came with the flow cytometer.  To continue analyzing the data she collected over the summer as part of her project she’ll need a different solution.

To give her a way to work with the FCS files I put together a quick R script that reads in the file, sets some event limits, and produces a nice plot.  With a little modification one could “gate” and count different regions.  The script uses the flowCore package to read in the FCS format files, and the hist2d command in gplots to make a reasonably informative plot.

library('flowCore')
library('gplots')

#### parameters ####

f.name <- 'file.name.goes.here'  # name of the file you want to analyze, file must have extension ".FCS"
sample.size <- 1e5               # number of events to plot, use "max" for all points
fsc.ll <- 1                      # FSC lower limit
ssc.ll <- 1                      # SSC lower limit
fl1.ll <- 1                      # FL1 lower limit (ex488/em536)

#### functions ####

## plotting function

plot.events <- function(fcm, x.param, y.param){
  hist2d(log10(fcm[,x.param]),
         log10(fcm[,y.param]),
         col = c('grey', colorRampPalette(c('white', 'lightgoldenrod1', 'darkgreen'))(100)),
         nbins = 200,
         bg = 'grey',
         ylab = paste0('log10(', y.param, ')'),
         xlab = paste0('log10(', x.param, ')'))
  
  box()
}

#### read in file ####

fcm <- read.FCS(paste0(f.name, '.FCS'))
fcm <- as.data.frame((exprs(fcm)))

#### analyze file and make plot ####

## eliminate values that are below or equal to thresholds you
## defined above

fcm$SSC[fcm$SSC <= ssc.ll|fcm$FSC <= fsc.ll|fcm$FL1 == fl1.ll] <- NA
fcm <- na.omit(fcm)

fcm.sample <- fcm

if(sample.size != 'max'){
  try({fcm.sample <- fcm[sample(length(fcm$SSC), sample.size),]},
      silent = T)
}

## plot events in a couple of different ways

plot.events(fcm, 'FSC', 'SSC')
plot.events(fcm, 'FSC', 'FL1')

## make a presentation quality figure

png(paste0(f.name, '_FSC', '_FL1', '.png'),
    width = 2000,
    height = 2000,
    pointsize = 50)

plot.events(fcm, 'FSC', 'FL1')

dev.off()

And here’s the final plot:

Posted in Computer tutorials, Research | Leave a comment

Microbial community segmentation with R

Posted on February 2, 2017 by Jeff

In my previous post I discussed our recent paper in ISME J, in which we used community structure and flow cytometry data to predict bacterial production.  The insinuation is that if you know community structure, and have the right measure of physiology, you can make a decent prediction of any microbial ecosystem function.  The challenge is that community structure data, which often has hundreds or thousands of dimensions (taxa, OTUs, etc.), is not easily used in straightforward statistical models.   Our workaround is to reduce the community structure data from many dimensions to a single categorical variable represented by a number.  We call this process segmentation.

You could carry out this dimension reduction with pretty much any clustering algorithm; you’re simply grouping samples with like community structure characteristics on the assumption that like communities will have similar ecosystem functions.  We  use the emergent self organizing map (ESOM), a neural network algorithm, because it allows new data to be classified into an existing ESOM.  For example, imagine that you are collecting a continuous time series of microbial community structure data.  You build an ESOM to segment your first few years of data, subsequent samples can be quickly classified into the existing model.  Thus the taxonomic structure, physiological, and ecological characteristics of the segments are stable over time.  There are other benefits to use an ESOM.  One is that with many samples (far more than we had in our study), the ESOM is capable of resolving patterns that many other clustering techniques cannot.

There are many ways to construct an ESOM.  I haven’t tried a Python-based approach, although I’m keen to explore those methods.  For the ISME J paper I used the Kohonen package in R, which has a nice publication that describes some applications and is otherwise reasonably well documented.  To follow this tutorial you can download our abundance table here.  Much of the inspiration, and some of the code for this analysis, follows the (retail) customer segmentation example given here.

For this tutorial you can download a table of the closest estimated genomes and closest completed genomes (analogous to an abundance table) here.  Assuming you’ve downloaded the data into your working directory, fire up Kohonen and build the ESOM.

## Kohonen needs a numeric matrix
edge.norm <- as.matrix(read.csv('community_structure.csv', row.names = 1))

## Load the library
library('kohonen')

## Define a grid.  The bigger the better, but you want many fewer units in the grid
## than samples.  1:5 is a good ballpark, here we are minimal.
som.grid <- somgrid(xdim = 5, ydim=5, topo="hexagonal")

## Now build the ESOM!  It is worth playing with the parameters, though in
## most cases you will want the circular neighborhood and toroidal map structure.
som.model.edges <- som(edge.norm, 
                 grid = som.grid, 
                 rlen = 100,
                 alpha = c(0.05,0.01),
                 keep.data = TRUE,
                 n.hood = "circular",
                 toroidal = T)

Congratulations!  You’ve just constructed your first ESOM.  Pretty easy.  You’ve effectively clustered the samples into the 25 units that define the ESOM.  You can visualize this as such:

plot(som.model.edges, type = 'mapping', pch = 19)

There are the 25 map units, with the toroid split and flattened into 2D.  Each point is a sample (row in the abundance table), positioned in the unit that best reflects its community structure.  I’m not going to go into any depth on the ESOM algorithm, which is quite elegant, but the version implemented in the Kohonen package is based on Euclidean distance.  How well each map unit represents the samples positioned within it is represented by the distance between the map unit and each sample.  This can be visualized with:

plot(som.model.edges, type = 'quality', pch = 19, palette.name = topo.colors)

Units with shorter distances in the plot above are better defined by the samples in those units than units with long distances.  What distance is good enough depends on your data and objectives.

The next piece is trickier because there’s a bit of an art to it.  At this point each sample has been assigned to one of the 25 units in the map.  In theory we could call each map unit a “segment” and stop here.  It’s beneficial however, to do an additional round of clustering on the map units themselves.  Particularly on large maps (which clearly this is not) this will highlight major structural features in the data.  Both k-means and hierarchical clustering work fairly well, anecdotally k-means seems to work better with smaller maps and hierarchical with larger maps, but you should evaluate for your data.  Here we’ll use k-means.  K-means requires that you specify the number of clusters in advance, which is always a fun chicken and egg problem.  To solve it we use the within-clusters sum of squares method:

wss.edges <- (nrow(som.model.edges$codes)-1)*sum(apply(som.model.edges$codes,2,var)) 
for (i in 2:15) {
  wss.edges[i] <- sum(kmeans(som.model.edges$codes, centers=i)$withinss)
}

plot(wss.edges,
     pch = 19,
     ylab = 'Within-clusters sum of squares',
     xlab = 'K')

Here’s where the art comes in.  Squint at the plot and try to decide the inflection point.  I’d call it 8, but you should experiment with whatever number you pick to see if it makes sense downstream.

We can make another plot of the map showing which map units belong to which clusters:

k <- 8
som.cluster.edges <- kmeans(som.model.edges$codes, centers = k)

plot(som.model.edges,
     main = '',
     type = "property",
     property = som.cluster.edges$cluster,
     palette.name = topo.colors)
add.cluster.boundaries(som.model.edges, som.cluster.edges$cluster)

Remember that the real shape of this map is a toroid and not a square.  The colors represent the final “community segmentation”; the samples belong to map units, and the units belong to clusters.  In our paper we termed these clusters “modes” to highlight the fact that there are real ecological properties associated with them, and that (unlike clusters) they support classification.  To get the mode of each sample we need to index the sample-unit assignments against the unit-cluster assignments.  It’s a little weird until you get your head wrapped around it:

som.cluster.edges$cluster[som.model.edges$unit.classif]
[1] 5 7 7 5 2 7 5 3 7 5 2 6 1 1 1 7 5 4 7 7 5 7 7 7 7 7 7 1 4 4 4 4 7 7 7 6 6 6 6 1 1 1 7 5 5 5 1 1 1 5 5 7 7 4 8 7 7 4 7 8
[61] 7 7 7 7 6 5 6 7 7 7 6 4 6 5 4 4 6 2 1 1 1 1 1 4 1 4 4 4

A really important thing to appreciate about these modes is that they are not ordered or continuous.  Mode 4 doesn’t necessarily have more in common with mode 5 say, than with mode 1.  For this reason it is important to treat the modes as factors in any downstream analysis (e.g. in linear modeling).  For our analysis I had a dataframe with bacterial production, chlorophyll concentration, and bacterial abundance, and predicted genomic parameters from paprica.  By saving the mode data as a new variable in the dataframe, and converting the dataframe to a zoo timeseries, it was possible to visualize the occurrence of modes, model the data, and test the pattern of modes for evidence of succession.  Happy segmenting!

Posted in Research | 6 Comments

New paper published in ISME Journal

Posted on January 20, 2017 by Jeff

I’m happy to report that a paper I wrote during my postdoc at the Lamont-Doherty Earth Observatory was published online today in the ISME Journal.  The paper, Bacterial community segmentation facilitates the prediction of ecosystem function along the coast of the western Antarctic Peninsula, uses a novel technique to “segment” the microbial community present in many different samples into a few groups (“modes”) that have specific functional, ecological, and genomic attributes.  The inspiration for this came when I stumbled across this blog entry on an approach used in marketing analytics.  Imagine that a retailer has a large pool of customers that it would like to pester with ads tailored to purchasing habits.  It’s too cumbersome to develop an individualized ad based on each customer’s habits, and it isn’t clear what combination of purchasing-habit parameters accurately describe meaningful customer groups.  Machine learning techniques, in this case emergent self-organizing maps (ESOMs), can be used to sort the customers in a way that optimizes their similarity and limits the risk of overtraining the model (including parameters that don’t improve the model).

In a 2D representation of an ESOM, the customers most like one another will be organized in geographically coherent regions of the map.  Hierarchical or k-means clustering can be superimposed on the map to clarify the boundaries between these regions, which in this case might represent customers that will respond similarly to a targeted ad.  But what’s really cool about this whole approach is that, unlike with NMDS or PCA or other multivariate techniques based on ordination, new customers can be efficiently classified into the existing groups.  There’s no need to rebuild the model unless a new type of customer comes along, and it is easy to identify when this occurs.

Back to microbial ecology.  Imagine that you have a lot of samples (in our case a five year time series), and that you’ve described community structure for these samples with 16S rRNA gene amplicon sequencing.  For each sample you have a table of OTUs, or in our case closest completed genomes and closest estimated genomes (CCGs and CEGs) determined with paprica.  You know that variations in community structure have a big impact on an ecosystem function (e.g. respiration, or nitrogen fixation), but how to test the correlation?  There are statistical methods in ecology that get at this, but they are often difficult to interpret.  What if community structure could be represented as a simple value suitable for regression models?

Enter microbial community segmentation.  Following the customer segmentation approach described above, the samples can be segmented into modes based on community structure with an Emergent Self Organizing Map and k-means clustering.  Here’s what this looks like in practice:

From Bowman et al. 2016.  Segmentation of samples based on bacterial community structure.  C-I show the relative abundance of CEGs and CCGs in each map unit.  This value was determined iteratively while the map was trained, and reflects the values for samples located in each unit (B).

This segmentation reduces the data for each sample from many dimensions (the number of CCG and CEG present in each samples) to 1.  This remaining dimension is a categorical variable with real ecological meaning that can be used in linear models.  For example, each mode has certain genomic characteristics:

From Bowman et al. 2016.  Genomic characteristics of modes (a and b), and metabolic pathways associated with taxa that account for most of the variations in composition between modes (d).

In panel a above we see that samples belonging to modes 5 and 7 (dominated by the CEG Rhodobacteraceae and CCG Dokdonia MED134, see Fig. 2 above) have the greatest average number of 16S rRNA gene copies.  Because this is a characteristic of fast growing, copiotrophic bacteria, we might also associate these modes with high levels of bacterial production.

Because the modes are categorical variables we can insert them right into linear models to predict ecosystem functions, such as bacterial production.  Combined with bacterial abundance and a measure of high vs. low nucleic acid bacteria, mode accounted for 76 % of the variance in bacterial production for our samples.  That’s a strong correlation for environmental data.  What this means in practice is; if you know the mode, and you have some flow cytometry data, you can make a pretty good estimate of carbon assimilation by the bacterial community.

For more on what you can do with modes (such as testing for community succession) check out the article!  I’ll post a tutorial on how to segment microbial community structure data into modes using R in a separate post.  It’s easier than you think…

Posted in Research | Leave a comment

paprica v0.4.0

Posted on January 8, 2017 by Jeff

I’m happy to announce the release of paprica v0.4.0.  This release adds a number of new features to our pipeline for evaluating microbial community and metabolic structure.  These include:

  • NCBI taxonomy information for each point of placement on the reference tree, including internal nodes.
  • Inclusion of the domain Eukarya.  This was a bit tricky and requires some further explanation.

The distribution of metabolic pathways, predicted during the creation of the paprica Eukarya database, across transcriptomes in the MMETSP.

Eukaryotic genomes are a totally different beast than their archaeal and bacterial counterparts.  First and foremost they are massive.  Because of these there aren’t very many completed eukaryotic genomes out there, particularly for singled celled eukaryotes.  While a single investigator can now sequence, assemble, and annotate a bacterial or archaeal genome in very little time, eukaryotic genomes still require major efforts by consortia and lots of $$.

One way to get around this scale problem is to focus on eukaryotic transcriptomes instead of genomes.  Because much of the eukaryotic genome is noncoding this greatly reduces sequencing volume.  Since there is no such thing as a contiguous transcriptome, this approach also implies that no assembly (beyond open reading frames) will be attempted.  The Moore Foundation-funded Marine Microbial Eukaryotic Transcriptome Sequencing Project (MMETSP) was an initial effort to use this approach to address the problem of unknown eukaryotic genetic diversity.  The MMETSP sequenced transcriptomes from several hundred different strains.  The taxonomic breadth of the strains sequenced is pretty good, even if (predictably) the taxonomic resolution is not.  Thus, as for archaea, the phylogenetic tree and metabolic inferences should be treated with caution.  For eukaryotes there are the additional caveats that 1) not all genes coded in a genome will be represented in the transcriptome 2) the database contains only strains from the marine environment and 3) eukaryotic 18S trees are kind of messy.  Considerable effort went into making a decent tree, but you’ve been warned.

Because the underlying data is in a different format, not all genome parameters are calculated for the eukaryotes.  18S gene copy number is not determined (and thus community and metabolic structure are not normalized), the phi parameter, GC content, etc. are also not calculated.  However, eukaryotic community structure is evaluated and metabolic structure inferred in the same way as for the domains bacteria and archaea:

./paprica-run.sh test.eukarya eukarya

As always you can install paprica v0.4.0 by following the instructions here, or you can use the virtual box or Amazon Web Service machine instance.

Posted in paprica | Leave a comment

paprica on the cloud

Posted on November 25, 2016 by Jeff

This is a quick post to announce that paprica, our pipeline to evaluate community structure and conduct metabolic inference, is now available on the cloud as an Amazon Machine Instance (AMI).  The AMI comes with all dependencies required to execute the paprica-run.sh script pre-installed.  If you want to use it for paprica-build.sh you’ll have to install pathway-tools and a few additional dependencies.  I’m new to the Amazon EC2 environment, so please let me know if you have any issues using the AMI.

If you are new to Amazon Web Services (AWS) the basic way this works is:

  • Sign up for Amazon EC2 using your normal Amazon log-in
  • From the AWS console, make sure that your region is N. Virginia (community AMI’s are only available in the region they were created in)
  • From your EC2 dashboard, scroll down to “Create Instance” and click “Launch Instance”
  • Now select the “Community AMIs”
  • Search for “paprica-ec2”, then select the AMI corresponding to the latest version of paprica (0.4.0 at the time of writing).
  • Choose the type of instance you would like to run the AMI on.  This is the real power of AWS; you can tailor the instance to the analysis you would like to run.  For testing choose the free t2.micro instance.  This is sufficient to execute the test files or run a small analysis (hundreds of reads).  To use paprica’s parallel features select an instance with the desired number of cores and sufficient memory.
  • Click “Review and Launch”, and finish setting up the instance as appropriate.
  • Log onto the instance, navigate to the paprica directory, execute the test file(s) as described in the paprica tutorial.  The AMI is not updated as often as paprica, so you may wish to reclone the github repository, or download the latest stable release.
Posted in paprica | 1 Comment

Antarctic Long Term Ecological Research

Posted on October 13, 2016 by Jeff
stuff

The Palmer and McMurdo Long Term Ecological Research (LTER) projects; separate but equal… at least in terms of interesting ecosystem dynamics if not in terms of biomass!

I’m very excited that our manuscript “Microbial community dynamics in two polar extremes: The lakes of the McMurdo Dry Valleys and the West Antarctic Peninsula Marine Ecosystem” has been published as an overview article in the journal BioScience.  The article belongs to a special issue comparing different ecological aspects of the two NSF-funded Long Term Ecological Research (LTER) sites in Antarctica.  I’m actually writing this post on my return trip from the first ever open science meeting of the International Long Term Ecological Research (ILTER) network at Kruger National Park in South Africa (an excellent place to ponder ecological questions).

This article had an odd genesis; the special issue was conceived by John Priscu, a PI with the McMurdo LTER project.  I was ensnared in the project along with Trista Vick-Majors, a graduate student with John Priscu (now a postdoctoral scholar at McGill University), shortly after starting my postdoc with Hugh Ducklow, PI on the Palmer LTER project.  The guidance we received was more or less “compare the McMurdo and Palmer LTERs”.  How exactly we should compare perennially ice-covered lakes in a polar desert to one of the richest marine ecosystems on the planet was left up to us.  Fortunately, microbial ecology lends itself to highly reductionist thinking.   This isn’t always helpful, but we reasoned that on a basal level the two ecosystems must function more or less the same.  Despite dramatically different physical settings, both environments host communities of phytoplankton (sometimes even similar taxonomic groups).  These convert solar energy into chemical energy and CO2 into organic carbon, thereby supporting communities of heterotrophic bacteria and grazers.

To look at the details of this we stretched the bounds of what constitutes an “overview article” and aggregated nearly two decades of primary production and bacterial production data collected by the McMurdo LTER, and over a decade of the same from the Palmer LTER.  By looking at the ratio of bacterial production to primary production we assessed how much carbon the heterotrophic bacterial community takes up relative to how much the phytoplankton community produces.

Some stuff

Figure from Bowman et al., 2016, BioScience.  A) Depth-integrated bacterial (BP) and primary production (PP) for the Palmer LTER study area and several lakes in Taylor Valley.  B)  The region occupied by the mean and standard deviation for discrete points (too many points to show).  C) The distribution of BP:PP for each site.

Typical marine values for this ratio are 1:10.  At a value of around 1:5 the carbon demands of heterotrophic bacteria are probably not met by phytoplankton production (the majority of carbon taken up by bacteria is lost through respiration and is not accounted for in the bacterial production assay).  Most of the lakes hover around 1:5, with values above this fairly common.  Lake Fryxell however, an odd lake at the foot of Canada Glacier, has values that often exceed 1:1!  Consistent with previous work on the lakes such high rates of bacterial production (relative to primary production) can only be met by a large external carbon subsidy.

Where does this external carbon come from?  Each summer the McMurdo Dry Valleys warm up enough that the various glaciers at the valley peripheries begin to melt.  This meltwater fuels chemoautotrophic bacterial communities where the glacier meets rock (the subglacial environment), and microbial mats in various streams and melt ponds.  Like microbial communities everywhere these bleed a certain amount of dissolved carbon (and particulate; DOC and POC) into the surrounding water.  Some of this carbon ends up in the lakes where it enhances bacterial production.

But external carbon subsidies aren’t the only part of the story.  Nutrients, namely phosphate and nitrate, are washed into the lakes as well.  During big melt years (such as the summer of 2001-2002 when a major positive SAM coupled to an El Nino caused unusually high temperatures) the lakes receives big pulses of relatively labile carbon but also inorganic nutrients and silt.  This odd combination has the effect of suppressing primary production in the near term through lowered light levels (all that silt), enhancing it in the long term (all those nutrients), and giving heterotrophic bacteria some high quality external carbon to feed on during the period that primary production is suppressed.  Or at least that’s how we read it.

Not a lake person?  How do things work over in the Palmer LTER?  One of the biggest ecological differences between Palmer and McMurdo is that the former has grazers (e.g. copepods, salps, and krill) and the latter does not, or at least not so many to speak off.  Thus an argument can be made that carbon dynamics at Palmer are driven (at least partially) by top-down controls (i.e. grazers), while at McMurdo they are dependent almost exclusively on bottom-up (i.e. chemical and physical) controls.

At times the difference between bacterial production and primary production is pretty extreme at Palmer.  In the summer of 2006 for example, bacterial production was only 3 % of primary production (see Fig. 4 in the publication), and the rate of primary production that summer was pretty high.  The krill population was also pretty high that year; at the top of their 4-year abundance cycle (see Saba et al. 2014, Nature Communications).  This is speculative, but I posit that bacterial production was low in part because a large amount of carbon was being transferred via krill to the higher trophic levels and away from bacteria.  This is a complicated scenario because krill can be good for bacteria; sloppy feeding produces DOC and krill excrete large amounts of reduced nitrogen and DOC.  Krill also build biomass and respire however, and their large fecal pellets sink quickly, these could be significant losses of carbon from the photic zone.

Antarctica is changing fast and in ways that are difficult to predict.  Sea ice seems to be growing in the east Antarctic as it is lost from the west Antarctic, and anomalous years buck this trend in both regions.  A major motivation for this special issue was to explore how the changing environment might drive ecological change.  I have to say that after spending a good portion of the (boreal) summer and fall thinking about this, some of that time from the vantage point of Palmer Station, I have no idea.  All of the McMurdo Lakes react differently to anomalous years, and Palmer as a region seems to react differently to each of abnormal year.  I think the krill story is an important concept to keep in mind here; ecological responses are like superimposed waveforms.  Picture a regularly occurring phenomenon like the El-Nino Southern Oscillation imposing a periodicity on sea ice cover, which we know has a strong influence on biology.  Add a few more oscillating waves from other physical processes.  Now start to add biological oscillations like the four-year krill abundance cycle.  Can we deconvolute this mess to find a signal?  Can we forecast it forward?  Certainly not with 10 years of data at one site and 20 years at the other (and we’re so proud of these efforts!).  Check back next century… if NSF funds these sites that long…

Many thanks to my co-authors for going the distance on this paper, particularly the lake people for many stimulating arguments.  I think limnology and oceanography are, conceptually, much less similar than lakes and oceans.

Posted in Research | 1 Comment

Creating a landmask for the West Antarctic Peninsula in R

Posted on September 4, 2016 by Jeff
presentation_graphic

Silicate concentration in µM at 1m depth during the 2014 Palmer LTER cruise.  This plot is lying to you.  The interpolations extend past islets and into landmasses.

This is going to be a pretty niche topic, but probably useful for someone out there.  Lately I’ve been working with a lot of geospatial data for the West Antarctic Peninsula.  One of the things that I needed to do was krig the data (krigging is a form of 2D interpolation, I’m using the pracma library for this).  Krigging is a problem near coastlines because it assumes a contiguous space to work in.  If there happens to be an island or other landmass in the way there is no way to represent the resulting discontinuity in whatever parameter you’re looking at.  Because of this I needed to find a way to mask the output.  This doesn’t really solve the problem, but at least it allows me to identify areas of concern (for example interpolation that extends across an isthmus, if there are sample points only on one side).

I’m krigging and building maps entirely inside R, which has somewhat immature packages for dealing with geospatial data.  The easiest masking solution would be to use filled polygons from any polygon format shapefile that accurately represents the coastline.  Unfortunately I couldn’t find an R package that does this correctly with the shapefiles that I have access too.  In addition, because of my downstream analysis it was better to mask the data itself, and not just block out landmasses in the graphical output.

Sharon Stammerjohn at the Institute of Arctic and Alpine Research pointed me to the excellent Bathymetry and Global Relief dataset produced by NOAA.  This wasn’t a complete solution to the problem but it got me moving in the right direction.  From the custom grid extractor at http://maps.ngdc.noaa.gov/viewers/wcs-client/ I selected a ETOPO1 (bedrock) grid along the WAP, with xyz as the output format.  If you’re struggling with shapefiles the xyz format is like a cool drink of water, being a 3-column matrix of longitude, latitude, and height (or depth).  For the purpose of creating the mask I considered landmass as any lat-long combination with height > 0.

There is one more really, really big twist to what I was trying to do, however.  The Palmer LTER uses a custom 1 km pixel grid instead of latitude-longitude.  It’s a little easier to conceptualize than lat-long given the large longitude distortions at high latitude (and the inconvenient regional convergence of lat-long values on similar negative numbers).  It is also a little more ecologically relevant, being organized parallel to the coastline instead of north to south.  Unfortunately this makes the grid completely incompatible with other Euclidean reference systems such as UTM.  So before I could use my xyz file to construct a land mask I needed to convert it to the line-station grid system used by the Palmer LTER.  If you’re working in lat-long space you can skip over this part.

grid

The Palmer LTER grid provides a convenient geospatial reference for the study area, but converting between line (y) and station (x) coordinates and latitude-longitude is non-trivial.

Many moons ago someone wrote a Matlab script to convert lat-long to line-station which you can find here.  Unfortunately I’m not a Matlab user, nor am I inclined to become one.  Fortunately it was pretty straightforward to copy-paste the code into R and fix the syntatic differences between the two languages.  Three functions in total are required:

## AUTHORS OF ORIGINAL MATLAB SCRIPT:
#   Richard A. Iannuzzi
#   Lamont-Doherty Earth Observatory
#   iannuzzi@ldeo.columbia.edu
#   based on: LTERGRID program written by Kirk Waters (NOAA Coastal Services Center), February 1997

## some functions that are used by the main function

SetStation <- function(e, n, CENTEREAST, CENTERNORTH, ANGLE){
  uu = e - CENTEREAST
  vv = n - CENTERNORTH
  z1 = cos(ANGLE)
  z2 = sin(ANGLE)
  NorthKm = (z1 * uu - z2 *vv) / 1000 + 600
  EastKm = (z2 * uu + z1 * vv) / 1000 + 40
  
  return(c(NorthKm, EastKm))
}

CentralMeridian <- function(iz){
  if(abs(iz) > 30){
    iutz = abs(iz) - 30
    cm = ((iutz * 6.0) -3.0) * -3600
  }
  else{
    iutz = 30 - abs(iz)
    cm = ((iutz * 6.0) +3.0) * +3600
  }
  return(cm)
}

GeoToUTM <- function(lat, lon, zone){

  axis = c(6378206.4,6378249.145,6377397.155,
          6378157.5,6378388.,6378135.,6377276.3452,
          6378145.,6378137.,6377563.396,6377304.063,
          6377341.89,6376896.0,6378155.0,6378160.,
          6378245.,6378270.,6378166.,6378150.)
  
  bxis = c(6356583.8,6356514.86955,6356078.96284,
          6356772.2,6356911.94613,6356750.519915,6356075.4133,
          6356759.769356,6356752.31414,6356256.91,6356103.039,
          6356036.143,6355834.8467,6356773.3205,6356774.719,
          6356863.0188,6356794.343479,6356784.283666,
          6356768.337303)
  
  ak0 = 0.9996
  radsec = 206264.8062470964  
  sphere = 9
  
  a = axis[sphere - 1]             # major axis size
  b = bxis[sphere - 1]             # minior axis size
  es = ((1-b^2/a^2)^(1/2))^2       # eccentricity squared
  slat = lat * 3600                # latitude in seconds
  slon = -lon * 3600               # longitude in seconds
  cm = 0                           # central meridian in sec
  iutz = 0
  
  cm = CentralMeridian(zone)       # call the function
  
  phi = slat/radsec
  dlam = -(slon - cm)/radsec
  epri = es/(1.0 - es)
  en = a/sqrt(1.0 - es * sin(phi)^2)
  t = tan(phi)^2
  c = epri * cos(phi)^2
  aa = dlam * cos(phi)
  s2 = sin(2.0 * phi)
  s4 = sin(4.0 * phi)
  s6 = sin(6.0 * phi)
  f1 = (1.0 - (es/4.)-(3.0*es*es/64.)-(5.0*es*es*es/256))
  f2 = ((3*es/8)+(3.0*es*es/32)+(45*es*es*es/1024))
  f3 = ((15*es*es/256)+(45*es*es*es/1024))
  f4 = (35*es*es*es/3072)
  em = a*(f1*phi - f2*s2 + f3*s4 - f4*s6)
  xx = ak0 * en * (aa + (1.-t+c) * aa^3/6 + (5 - 18*t + t*t + 72*c-58*epri)* aa^5/120) + 5e5
  yy = ak0 * (em + en * tan(phi) *((aa*aa/2) + (5-t+9*c+4*c*c)*aa^4/24 + (61-58*t +t*t +600*c - 330*epri)* aa^6/720))
  
  if(zone < 0 | slat < 0){
    yy = yy + 1e7
  }
  
  return(c(xx, yy))
}

## This function actually works with your data

ll2gridLter <- function(inlat, inlon){
  NorthKm = 0           # initialize
  EastKm = 0            # initialize
  zone = -20            # set zone (for LTER region, I think)
  ANGLE = -50 * pi / 180
  CENTEREAST = 433820.404        # eastings for station 600.040
  CENTERNORTH = 2798242.817      # northings for station 600.040
  
  # take latitude longitude and get station
  x.y = GeoToUTM(inlat, inlon, zone)
  NorthKm.EastKm = SetStation(x.y[1], x.y[2], CENTEREAST, CENTERNORTH, ANGLE)
  return(NorthKm.EastKm)
}

Once the functions are defined I used them to convert the lat/long coordinates in the xyz file to line-station.

## Read in xyz file.
lat.long.depth <- read.table('etopo1_bedrock.xyz', header = F, col.names = c('long', 'lat', 'depth'))

## Limit to points above sea level.
lat.long.land <- lat.long.depth[which(lat.long.depth$depth >= 0),]

## Create a matrix to hold the output.
line.station.land <- matrix(ncol = 3, nrow = length(lat.long.land$long))
colnames(line.station.depth) <- c('line', 'station', 'depth')

## Execute the ll2gridLter function on each point.  Yes, I'm using a loop to do this.
for(i in 1:length(lat.long.land$long)){
  line.station.land[i,] <- c(ll2gridLter(lat.long.land$lat[i], lat.long.land$long[i]), lat.long.land$depth[i])
  print(paste(c(i, line.station.land[i,])))
}

## Write out the matrix.
write.csv(line.station.land, 'palmer_grid_landmask.csv', row.names = F, quote = F)

At this point I had a nice csv file with line, station, and elevation.  I was able to read this into my existing krigging script and convert into a mask.

## Read in csv file.
landmask <- read.csv('palmer_grid_landmask.csv')

## Limit to the lines and stations that I'm interested in.
landmask <- landmask[which(landmask[,1] <= 825 & landmask[,1] >= -125),]
landmask <- landmask[which(landmask[,2] <= 285 & landmask[,2] >= -25),]

## Interpolated data is at 1 km resolution, need to round off
## to same resolution here.
landmask.expand <- cbind(ceiling(landmask[,1]), ceiling(landmask[,2]))

## Unfortunately this doesn't adequately mask the land.  Need to expand the size of each
## masked pixel 1 km in each direction.
landmask.expand <- rbind(landmask.expand, cbind(floor(landmask[,1]), floor(landmask[,2])))
landmask.expand <- rbind(landmask.expand, cbind(ceiling(landmask[,1]), floor(landmask[,2])))
landmask.expand <- rbind(landmask.expand, cbind(floor(landmask[,1]), ceiling(landmask[,2])))
landmask.expand <- unique(landmask.expand)

I’m not going to cover how I did the krigging in this post.  My krigged data is in matrix called temp.pred.matrix with colnames given by ‘x’ followed by ‘station’, as in x20 for station 20, and row names ‘y’ followed by ‘line’, as in y100 for line 100.  To convert interpolated points that are actually land to NA values I simply added this line to my code:

temp.pred.matrix[cbind(paste0('y', landmask.expand[,1]), paste0('x', landmask.expand[,2] * -1))]

Here’s what the krigged silicate data looks like after masking.

stuff

Silicate concentration in µM at 1m depth during the 2014 Palmer LTER cruise after masking the underlying data.

Excellent.  The masked area corresponds with known landmasses; that’s Adelaide Island (home of Rothera Station) coming in at the eastern portion of Line 300, and various islets and the western edge of the Antarctic Peninsula to the northeast.  At this point erroneous data has been eliminated from the matrix.  Annual inventories of silicate and other nutrients can be more accurately calculated form the data and our eyes are not drawn to interesting features in the interpolation that have no chance of reflecting reality because they are over land.  The white doesn’t look that appealing to me in this plot however, so I masked the land with black by adding points to the plot.  Again, I’m not going to show the whole plotting routine because some variables would require a lengthy explanation about how my larger dataset is structured.  The plot was created using imagep in the oce package.  One thing to be aware of is that this command automatically transposes the matrix, thus mask needs to be transposed as well.  I did this manually in the following line:

## Show masked points as black squares.
points(landmask.expand[,1] ~ {-1 * landmask.expand[,2]}, pch = 15, cex = 0.6)

And the final plot:

si

Silicate concentration in µM at 1m depth during the 2014 Palmer LTER cruise after masking the underlying data and adding points to indicate the masked area.

 

 

Posted in Research | Leave a comment

Astrobiology Primer v2

Posted on August 20, 2016 by Jeff

The long-awaited version 2 of the Astrobiology Primer was published (open access) yesterday in the journal Astrobiology.  I’m not sure who first conceived of the Astrobiology Primer, first published in 2006, but the v2 effort was headed by co-lead editors Shawn Domagal-Goldman at NASA and Katherine Wright  at the UK Space Agency.  The Primer v2 was a long time in coming; initial section text was submitted back in early 2011!  The longer these projects go on, the easy it is for them to die.  Many thanks to Shawn and Katherine for making sure that this didn’t happen!

The downside of course, is that the primer ran the risk of being outdated before it even went into review.  This was mitigated somewhat by the review process itself, and authors did have a chance to update their various sections.  Some sections are more stable than others; the section that I wrote with Shawn McGlynn (now at the Tokyo Institute of Technology) on how life uses energy for example, covers some fairly fundamental ground and is likely to stand the test of time.  Less so for sections that cover planetary processes in and outside of our solar system; paradigms are being broken in planetary science as fast as they form!

The Astrobiology Primer is a very valuable document because it takes a complicated and interdisciplinary field of study and attempts to summarize it for a broad audience.  Most of the Primer should be accessible to anyone with a basic grasp of science.  I wonder if it could even serve as a model for other disciplines.  What if the junior scientists in every discipline (perhaps roughly defined by individual NSF or NASA programs) got together once every five years to write an open-access summary of the major findings in their field?  This might provide a rich and colorful counterpoint to the valuable but often [dry, esoteric, top-down?  There’s an adjective that I’m searching for here but it escapes me] reports produced by the National Academies.

The co-lead editors were critical to the success of the Primer v2.  I haven’t explicitly asked Shawn and Kaitlin if they were compensated in any way for this activity – perhaps they rolled some of this work into various fellowships and such over the years.  More likely this was one more extracurricular activity carried out on the side.  Such is the way science works, and the lines are sufficiently blurred between curricular and extracurricular that most of us don’t even look for them anymore.  In recognition of this, and to speed the publication and heighten the quality of a future Primer v3, it would be nice to see NASA produce a specific funding call for a (small!) editorial team.  Three years of partial salary and funding for a series of writing workshops would make a huge difference in scope and quality.

Posted in Research | Leave a comment

How I learned to stop worrying and love subsampling (rarifying)

Posted on July 12, 2016 by Jeff

I have had the 2014 paper “Waste Not, Want Not: Why Rarefying Microbiome Data is Inadmissable” by McMurdie and Holmes sitting on my desk for a while now.  Yesterday I finally got around to reading it and was immediately a little skeptical as a result of the hyperbole with which they criticized the common practice of subsampling* libraries of 16S rRNA gene reads during microbial community analysis.  The logic of that practice proceeds like this:

To develop 16S rRNA gene data (or any other marker gene data) describing a microbial community we generally collect an environmental sample, extract the DNA, amplify it, and sequence the amplified material.  The extract might contain the DNA from 1 billion microbial cells present in the environment.  Only a tiny fraction (<< 1 %) of these DNA molecules actually get sequenced; a “good” sequence run might contain only tens of thousands of sequences per sample.  After quality control and normalizing for multiple copies of the 16S rRNA gene it will contain far fewer.  Thus the final dataset contains a very small random sample of the original population of DNA molecules.

Most sequence-based microbial ecology studies involve some kind of comparison between samples or experimental treatments.  It makes no sense to compare the abundance of taxa between two datasets of different sizes, as the dissimilarity between the datasets will appear much greater than it actually is.  One solution is to normalize by dividing the abundance of each taxa by the total reads in the dataset to get their relative abundance.  In theory this works great, but has the disadvantage that it does not take into account that a larger dataset has sampled the original population of DNA molecules deeper.  Thus more rare taxa might be represented.  A common practice is to reduce the amount of information present in the larger dataset by subsampling to the size of the smaller dataset.  This attempts to approximate the case where both datasets undertake the same level of random sampling.

McMurdie and Holmes argue that this approach is indefensible for two common types of analysis; identifying differences in community structure between multiple samples or treatments, and identifying differences in abundance for specific taxa between samples and treatments.  I think the authors do make a reasonable argument for the latter analysis; however, at worst the use of subsampling and/or normalization simply reduces the sensitivity of the analysis.  I suspect that dissimilarity calculations between treatments or samples using realistic datasets are much less sensitive to reasonable subsampling than the authors suggest.  I confess that I am (clearly) very far from being a statistician and there is a lot in McMurdie and Holmes, 2014 that I’m still trying to digest.  I hope that our colleagues in statistics and applied mathematics continue optimizing these (and other) methods so that microbial ecology can improve as a quantitative science.  There’s no need however, to get everyone spun up without a good reason.  To try and understand if there is a good reason, at least with respect to my data and analysis goals, I undertook some further exploration.  I would strongly welcome any suggestions, observations, and criticisms to this post!

My read of McMurdi and Homes is that the authors object to subsampling because it disregards data that is present that could be used by more sophisticated (i.e. parametric) methods to estimate the true abundance of taxa.  This is true; data discarded from larger datasets does have information that can be used to estimate the true abundance of taxa among samples.  The question is how much of a difference does it really make?  McMurdi and Holmes advocate using methods adopted from transcriptome analysis.  These methods are necessary for transcriptomes because 1) the primary objective of the study is not usually to see if one transcriptome is different from another, but which genes are differentially expressed and 2) I posit that the abundance distribution of transcript data is different than the taxa abundance distribution.  An example of this can be seen in the plots below.

Taxon abundance, taken from a sample that I’m currently working with.

Transcript abundance by PFAM, taken from a sample (selected at random) from the MMETSP.

In both the cases the data is log distributed, with a few very abundant taxa or PFAMs and many rare ones.  What constitutes “abundant” in the transcript dataset however, is very different than “abundant” in the community structure dataset.  The transcript dataset is roughly half the size (n = 9,000 vs. n = 21,000), nonetheless the most abundant transcript has an abundance of 209.  The most abundant OTU has an abundance of 6,589, and there are several very abundant OTUs.  Intuitively this suggests to me that the structure of the taxon dataset is much more reproducible via subsampling than the structure of the transcript dataset, as the most abundant OTUs have a high probability of being sampled.  The longer tail of the transcript data contributes to this as well, though of course this tail is controlled to a large extent by the classification scheme used (here PFAMs).

To get an idea of how reproducible the underlying structure was for the community structure and transcript data I repeatedly subsampled both (with replacement) to 3000 and 1300 observations, respectively.  For the community structure data this is about the lower end of the sample size I would use in an actual analysis – amplicon libraries this small are probably technically biased and should be excluded (McMurdi and Holmes avoid this point at several spots in the paper).  The results of subsampling are shown in these heatmaps, where each row is a new subsample.  For brevity only columns summing to an abundance > 100 are shown.

otu_abundance_heat trans_abundance_heatIn these figures the warmer colors indicate a higher number of observations for that OTU/PFAM.  The transcript data is a lot less reproducible than the community structure data; the Bray-Curtis dissimilarity across all iterations maxes out at 0.11 for community structure and 0.56 for the transcripts.  The extreme case would be if the data were normally distributed (i.e. few abundant and few rare observations, many intermediate observations).  Here’s what subsampling does to normally distributed data (n = 21,000, mean = 1000, sd = 200):

otu_random_abundance

If you have normally distributed data don’t subsample!

For the rest of us it seems that for the test dataset used here community structure is at least somewhat reproducible via subsampling.  There are differences between iterations however, what does this mean in the context of the larger study?

The sample was drawn from a multiyear time series from a coastal marine site.  The next most similar sample (in terms of community composition and ecology) was, not surprisingly, a sample taken one week later.  By treating this sample in an identical fashion, then combining the two datasets, it was possible to evaluate how easy it is to tell the two samples apart after subsampling.  In this heatmap the row colors red and black indicate iterations belonging to the two samples:

Clustering of repeated subsamplings from two similar samples.  Sample identity is given by the red or black color along the y-axis.

 

As this heatmap shows, for these two samples there is perfect fidelity.  Presumably with very similar samples this would start to break down, determining how similar samples need to be before they cannot be distinguished at a given level of subsampling would be a useful exercise.  The authors attempt to do this in Simlation A/Figure 5 in the paper, but it isn’t clear to me why their results are so poor – particularly given very different sample types and a more sophisticated clustering method than I’ve applied here.

As a solution – necessary for transcript data, normally distributed data, or for analyses of differential abundance, probably less essential for comparisons of community structure – the authors propose a mixture model approach that takes in account variance across replicates to estimate “real” OTU abundance.  Three R packages that can do this are mentioned; edgeR, DESeq, and metagenomeSeq.  The problem with these methods – as I understand them – is that they require experimental replicates.  According to the DESeq authors, technical replicates should be summed, and samples should be assigned to treatment pools (e.g. control, treatment 1, treatment 2…).  Variance is calculated within each pool and this is used to to model differences between pools.  This is great for a factor-based analysis, as is common in transcriptome analysis or human microbial ecology studies.  If you want to find a rare, potentially disease-causing strain differently present between a healthy control group and a symptomatic experimental group for example, this is a great way to go about it.

There are many environmental studies for which these techniques are not useful however, as it may be impractical to collect experimental replicates.  For example it is both undesirable and impractical to conduct triplicate or even duplicate sampling in studies focused on high-resolution spatial or temporal sampling.  Sequencing might be cheap now, but time and reagents are not.  Some of the methods I’m working with are designed to aggregate samples into higher-level groups – at which point these methods could be applied by treating within-group samples as “replicates” – but this is only useful if we are interested in testing differential abundance between groups (and doesn’t solve the problem of needing to get samples grouped in the first place).

These methods can be used to explore differential abundance in non-replicated samples, however, they are grossly underpowered when used without replication.  Here’s an analysis of differential abundance between the sample in the first heatmap above and its least similar companion from the same (ongoing) study using DESeq.  You can follow along with any standard abundance table where the rows are samples and the columns are variables.

library(DESeq)

## DESeq wants to oriented opposite how community abundance
## data is normally presented (e.g. to vegan)
data.dsq <- t(data)

## DESeq requires a set of conditions which are factors.  Ideally
## this would be control and treatment groups, or experimental pools
## or some such, but we don't have that here.  So the conditions are
## unique column names (which happen to be dates).
conditions <- factor(as.character(colnames(data.dsq)))

## As a result of 16S rRNA gene copy number normalization abundance
## data is floating point numbers, convert to integers.
data.dsq <- ceiling(data.dsq, 0)

## Now start working with DESeq.
data.ct <- newCountDataSet(as.data.frame(data.dsq), conditions = conditions)
data.size <- estimateSizeFactors(data.ct)

## This method and sharing mode is required for unreplicated samples.
data.disp <- estimateDispersions(data.size, method = 'blind', sharingMode="fit-only")

## And now we can execute a test of differential abundance for the
## two samples used in the above example.
test <- nbinomTest(data.disp, '2014-01-01', '2014-01-06')
test <- na.omit(test)

## Plot the results.
plot(test$baseMeanA, test$baseMeanB,
     #log = 'xy',
     pch = 19,
     ylim = c(0, max(test$baseMeanB)),
     xlim = c(0, max(test$baseMeanA)),
     xlab = '2014-01-01',
     ylab = '2009-12-14')

abline(0, 1)

## If we had significant differences we could then plot them like
## this:
points(test$baseMeanA[which(test$pval < 0.05)],
       test$baseMeanB[which(test$pval < 0.05)],
       pch = 19,
       col = 'red')

binomial

 

As we would expect there are quite a few taxa present in high abundance in one sample and not the other, however, none of the associated p-values are anywhere near significant.  I’m tempted to try to use subsampling to create replicates, which would allow an estimate of variance across subsamples and access to greater statistical power.  This is clearly not as good as real biological replication, but we have to work within the constraints of our experiments, time, and funding…

*You might notice that I’ve deliberately avoided using the terms “microbiome” and “rarefying” here.  In one of his comics Randall Munroe asserted that the number of made up words in a book is inversely proportional to the quality of the book, similarly I strongly suspect that the validity of a sub-discipline is inversely proportional to the number of jargonistic words that community makes up to describe its work.  As a member of said community, what’s wrong with subsampling and microbial community??

Posted in Research | 2 Comments