Ecuador update

Today we took the last sample of our Ecuador field effort, though we have a few days left in-country. Right now we are in the town of Mompiche, just down the coast from our second field site near Muisne. Tomorrow we’ll be sorting out gear and getting ready for a few days of meetings in Guayaquil. Then its time to fly home and start working up some data! I’m too tired to write a coherent post on the last two (intensive!) weeks of sampling, but here’s a photo summary of our work in the Cayapas-Mataje Ecological Reserve, where we spent the bulk of our time. Enjoy!

Pelicans congregating along a front in the Cayapas-Mataje estuary.

The town of Valdez, near the mouth of the Cayapas-Mataje estuary.  The Reserve is right on the border with Columbia, and up until a few years ago Valdez had a reputation as a trafficking hub.  Drug trafficking is still a problem throughout the Reserve, but with the conflict with FARC more or less over I understand that tension in the local communities has gone way down.  Valdez seems okay, and the people we met there were friendly.

Another view of Valdez.

Shrimp farm in the Cayapas-Mataje Ecological Reserve.  You can’t build a new farm in the Reserve, but old shrimp farms were grandfathered in when the Reserve was created.

The R/V Luz Yamira, at its home port of Tambillo.  Tambillo was a vibrant, friendly little town where we spent a bit of time.  The town is struggling to hold onto its subsistence fishing lifestyle in the face of declining fish stocks.

ADCP, easy as 1-2-3

Birds of a feather…

Morning commute from the city of San Lorenzo.

The Cayapas-Mataje Ecological Reserve has the tallest mangrove trees in the world.

I took this picture just for the good people at UCSD risk management.

Team Cayapas-Mataje.  From left; Jessie, Natalia, and Jeff.  We are standing in front of Jessie’s house in Tambillo.  Many thanks to Jessie and his wife for letting us stay a night and get away from the craziness of San Lorenzo!

A very full car ready to head to Muisne.  Its a good thing Natalia and I are both fairly short.

A large shrimp farm near Muisne.  The area around Muisne has been almost entirely deforested for shrimp aquaculture.  By comparing this area with the more pristine Cayapas-Mataje, we hope to better understand the biogeochemical consequences of coastal land-use change.

Posted in Uncategorized | Leave a comment

Initial field effort in Ecuador

As any reader of this blog will know, most of the research in the Bowman Lab is focused on polar microbial ecology. Although focusing a research program on a set of geographically-linked environments does have advantages, primarily the ability to spend more time thinking in depth about them, there is I think something lost with this approach. Insights are often gained by bringing a fresh perspective to a new study area, or applying lessons learned in such an area to places that one has studied for years. With this in mind lab member Natalia Erazo and I are launching a new field effort in coastal Ecuador. Natalia is an Ecuadorean native, and gets credit for developing the idea and sorting out the considerable logistics for this initial effort. Our first trip is very much a scouting effort, but will carry out some sampling in the Cayapas-Mataje Ecological Reserve and near the town of Muisne. Depending on funding we hope to return during the rainy season in January-February for a more intensive effort.

Dense coastal forest in the Cayapas-Mataje Ecological Reserve.  Photo Natalia Erazo.

Our primary objective is to understand the role of mangrove forests in coastal biogeochemical cycling. Mangroves are salt-tolerant trees that grow in tropical and sub-tropical coastal areas around the world. They are known to provide a range of positive ecosystem functions; serving as fish habitat, stabilizing shorelines, and providing carbon and nutrient subsidies to the coastal marine environment. Globally mangroves are under threat. The population density of many tropical coastal areas is increasing, and that inevitably leads to land-use changes (such as deforestation) and a loss of economic services – the social and economic benefits of an ecosystem – as economic activity ramps up. The trick to long-term sustainability is to maintain ecosystem services during economic development, allowing standards of living to increase in the short term without a long-term economic loss resulting from ecological failure (such as the collapse of a fishery or catastrophic coastal flooding). This is not easily done, and requires a much better understanding of what functions exactly, specific at-risk components of an ecosystem provide than we often have.

One particular land-use threat to the mangrove ecosystem is shrimp aquaculture. Mangrove forests in Ecuador and in other parts of the world have been deforested on a massive scale to make room for shrimp aquaculture ponds. In addition to scaling back any ecosystem functions provided by the mangrove forest, the shrimp aquaculture facilities are a source of nutrients in the form of excrement and excess feed. On this trip we will try to locate estuaries more and less perturbed by aquaculture. By comparing nutrient and carbon concentrations, sediment load, and microbial community structure between these areas, we will gain a preliminary understanding of what happens to the coastal ecosystem when mangroves are removed and aquaculture facilities are installed in their place.

Our first stop on this search will be San Lorenzo, a small city in the Cayapas-Mataje Ecological Reserve near the border with Columbia. I’m extremely excited to visit the Reserve, which has the distinction of hosting the tallest mangrove trees anywhere on Earth. We may some meager internet access in San Lorenzo, so I’ll try to update this blog as I’m able. Because of the remote nature of some of our proposed field sites we’ll have a Garmin InReach satellite messenger with us. We plan to leave the device on during our field outings, you can track our location in real time on the map below! The Garmin map interface is a bit cludgey; you should ignore the other Scripps Institution of Oceanography users listed on the side panel as I can’t seem to make them disappear.

Posted in Uncategorized | Leave a comment

So you want to use your computer for science…

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

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

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

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

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

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
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

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

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 “” and type:


echo "hello world!"
# hey, this line doesn't do anything!
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

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$ ./
hello world!

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

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$
-bash: command not found

The command would execute if you typed temp/, 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

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:

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$
hello world!

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

Microbial session at POLAR 2018 in Davos

Davos, Switzerland, site of the POLAR2018 conference.  Image from

With colleagues Maria Corsaro, Eric Collins, Maria Tutino, Jody Deming, and Julie Dinasquet I’m convening a session on polar microbial ecology and evolution at the upcoming POLAR2018 conference in Davos, Switzerland.  Polar2018 is shaping up to be a unique and excellent conference; for the first time (I think) the major international Arctic science organization (IASC) is joining forces with the major international Antarctic science organization (SCAR) for a joint meeting.  Because many polar scientists specialize in one region, and thus have few opportunities to learn from the other, a joint meeting will be a great opportunity to share ideas.

Here’s the full-text description for our session.  If your work is related to microbial evolution, adaption, or ecological function in the polar regions I hope you’ll consider applying!

Posted in Uncategorized | Leave a comment

Denitrifying bacteria and oysters

The eastern oyster.  Denitrification never tasted so good!  Photo from

I’m happy to be co-author on a study that was just published by Ann Arfken, a PhD student at the Virginia Institute for Marine Science (VIMS).  The study evaluated the composition of the microbial community associated with eastern oyster Crassostrea virginica to determine if oysters play a role in denitrification.  Denitrification is an ecologically significant anaerobic microbial metabolism.  In the absence of oxygen certain microbes use other oxidized compounds as electron acceptors.  Nitrate (NO3) is a good alternate electron acceptor, and the process of reducing nitrate to nitrite (NO2), and ultimately to elemental nitrogen (N2), is called denitrification.  Unfortunately nitrate is needed by photosynthetic organisms like plants and phytoplankton, so the removal of nitrate can be detrimental to ecosystem health.  Oxygen is easily depleted in the guts of higher organisms by high rates of metabolic activity, creating a niche for denitrification and other anaerobic processes.

Predicted relative abundance (paprica) as a function of measured (qPCR) relative abundance of nosZI genes.  From Arfken et al. 2017.

To evaluate denitrification in C. virginica, Arfken et al. coupled actual measurements of denitrification in sediments and oysters with an analysis of microbial community structure in oyster shells and digestive glands.  We then used paprica with a customized database to predict the presence of denitrification genes, and validated the predictions with qPCR.

I was particularly happy to see that the qPCR results agreed well with the paprica predictions for the nosZ gene, which codes for the enzyme responsible for reducing nitrous oxide (N2O) to N2.  I believe this is the first example of qPCR being used to validate metabolic inference – which currently lacks a good method for validation.  Surprisingly however (at least to me), denitrification in C. virginica was largely associated with the oyster shell rather than the digestive gland.  We don’t really know why this is.  Arfken et al. suggests rapid colonization of the oyster shell by denitrifying bacteria that also produce antibiotic compounds to exclude predators, but further work is needed to demonstrate this!

Posted in Uncategorized | 3 Comments

Analyzing flow cytometry data with R

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.


#### parameters #### <- ''  # 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 <- function(fcm, x.param, 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, ')'))

#### read in file ####

fcm <- read.FCS(paste0(, '.FCS'))
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, 'FSC', 'SSC'), 'FSC', 'FL1')

## make a presentation quality figure

png(paste0(, '_FSC', '_FL1', '.png'),
    width = 2000,
    height = 2000,
    pointsize = 50), 'FSC', 'FL1')

And here’s the final plot:

Posted in Research | Leave a comment

OPAG comes to Scripps Institution of Oceanography

Artist’s rendition of the Dawn spacecraft approaching Ceres.  Original image can be found here.

I’m excited to be hosting the fall meeting of NASA’s Outer Planets Assessment Group (OPAG) here at Scripps Institution of Oceanography in September.  For planetary scientists at UCSD, SDSU, USD, and other institutions in the greater San Diego area, if you’ve never attended the OPAG meeting here’s your chance!  The meeting will be September 6 and 7 at the Samual H. Scripps Auditorium.  More details can be found here.

What are assessment groups, and what is OPAG specifically?  The assessment groups are an excellent NASA innovation to encourage dialogue within the scientific community, and between the scientific community and NASA HQ.  There’s usually a little tense dialogue – in a good way – between these two ends of the scientific spectrum.  I often wish NSF had a similar open-format dialogue with its user community!  The form of the OPAG meeting is 15 or 30 minute presentations on a range of topics relevant to the community.  These are often mission updates, planning updates for future missions, and preliminary results from the analysis of mission data.  NASA has quite a few assessment groups, ranging from the Small Body Assessment Group (SBAG – probably the AG with the catchiest acronym) to the Mars Assessment Group (MPAG).  OPAG covers everything in the solar system further from the sun than Mars.  If that covers your favorite planetary body, come and check it out!

It’s traditional to have a public evening lecture with the OPAG meeting.  For the upcoming meeting the lecture will be given at 7 pm on September 6 at the Samuel Scripps Auditorium by my friend and colleague Britney Schmidt, a planetary scientist from Georgia Tech and an expert on Europa and on Antarctic ice sheets.  Why and how one can develop that dual expertise will certainly be made clear in her talk.  There is no cost to attend, but an RSVP is recommended.  You can find more details and RSVP here.

Posted in Uncategorized | Leave a comment

Analyzing “broad” data with glmnet

Very often environmental datasets contain far fewer observations than we would like, and far more variables that might influence these observations.  This is the case for one of my projects in the Palmer LTER: I have annual observations of an environmental indicator (n = 19) and far more than 19 environmental conditions (predictors, p) that might influence the indicator.  This is a broad dataset (n >> p), and the problem can’t be solved with classic multiple regression (for which the number of observations must exceed the number of predictors).  Enter the R package glmnet.  Glmnet is an implementation of lasso, ridge, and elastic-net regression.  There are a limited number of glmnet tutorials out there, including this one, but I couldn’t find one that really provided a practical start to end guide.  Here’s the method that I came up with for using glmnet to select variables for multiple regression.  I’m not an expert in glmnet, so don’t take this as a definitive guide.  Comments are welcome!

First, let’s get some data.  The data are various climate indices with a 13- month time lag as well as sea ice extent, open water extent, and fractional open water extent.  There are a couple of years missing from this dataset, and we need to remove one more year for reasons that aren’t important here.

env.params <- read.csv('', header = T, row.names = 1)
env.params <- env.params[row.names(env.params) != '1994',]

To keep things simple we’ll create our response variable by hand. Nevermind what the response variable is for now, you can read our paper when it comes out 🙂

response <- c(0.012, 0.076, 0.074, 0.108, 0.113, 0.154, 0.136, 0.183, 0.210, 0.043, 0.082, 0.092, 0.310, 0.185, 0.436, 0.357, 0.472, 0.631, 0.502)

One thing to note at this point is that, because we have so few observations, we can’t really withhold any to cross-validate the glmnet regression. That’s a shame because it means that we can’t easily optimize the alpha parameter (which determines whether glmnet uses lasso, elastic net, or ridge) as was done here. Instead we’ll use alpha = 0.9. This is to make the analysis a little bit elastic-netty, but still interpretable.

In my real analysis I was using glmnet to identify predictors for lots of response variables.  It was therefor useful to write a function to generalize the several commands needed for glmnet.  Here’s the function:


## The function requires a matrix of possible predictors, a vector with the response variable,
## the GLM family used for the model (e.g. 'gaussian'), the alpha parameter, and "type.measure".
## See the documentation on cv.glmnet for options.  Probably you want "deviance". 

get.glmnet <- function(predictors, response, family, alpha, type.measure){
  glmnet.out <- glmnet(predictors, response, family = family, alpha = alpha) <- cv.glmnet(predictors, response, family = family, alpha = alpha, type.measure = type.measure)
  ## Need to find the local minima for lambda.
  lambda.lim <-$lambda[which.min($cvm)]
  ## Now identify the coefficients that correspond to the lambda minimum.
  temp.coefficients.i <- coefficients(glmnet.out, s = lambda.lim)@i + 1 # +1 to account for intercept
  ## And the parameter names...
  temp.coefficients.names <- coefficients(glmnet.out, s = lambda.lim)@Dimnames[[1]][temp.coefficients.i]
  temp.coefficients <- coefficients(glmnet.out, s = lambda.lim)@x
  ## Package for export.
  temp.coefficients <- rbind(temp.coefficients.names, temp.coefficients)

Phew!  Okay, let’s try to do something with this.  Note that glmnet requires a matrix as input, not a dataframe.

response.predictors <- get.glmnet(data.matrix(env.params), response, 'gaussian', alpha, 'deviance')
      temp.coefficients.names temp.coefficients      
 [1,] "(Intercept)"           "0.496410195525042"    
 [2,] "ao.Apr"                "0.009282064516813"    
 [3,] "ao.current"            "-0.0214919836174853"  
 [4,] "pdo.Mar"               "-0.0568728879266135"  
 [5,] "pdo.Aug"               "-0.00881845994191182" 
 [6,] "pdo.Dec"               "-0.0321738320415234"  
 [7,] "ow.Dec"                "1.92231198892721e-06" 
 [8,] "ow.frac.current"       "0.207945851122607"    
 [9,] "ice.Sep"               "-2.29621552264475e-06"

So glmnet thinks there are 9 predictors.  Possible, but I’m a little suspicious.  I’d like to see this in the more familiar GLM format so that I can wrap my head around the significance of the variables.  Let’s start by building a null model of all the predictors.

null.model <- lm(response ~ ao.Apr + ao.current + pdo.Mar + pdo.Aug + pdo.Dec + ow.Dec + ow.frac.current + ice.Sep, data = env.params)
> summary(null.model)

lm(formula = response ~ ao.Apr + ao.current + pdo.Mar + pdo.Aug + 
    pdo.Dec + ow.Dec + ow.frac.current + ice.Sep, data = env.params)

     Min       1Q   Median       3Q      Max 
-0.14304 -0.03352 -0.01679  0.05553  0.13497 

                  Estimate Std. Error t value Pr(>|t|)  
(Intercept)      2.052e+00  1.813e+00   1.132   0.2841  
ao.Apr           2.226e-02  3.702e-02   0.601   0.5610  
ao.current      -3.643e-02  1.782e-02  -2.044   0.0681 .
pdo.Mar         -7.200e-02  3.215e-02  -2.240   0.0490 *
pdo.Aug         -1.244e-02  2.948e-02  -0.422   0.6821  
pdo.Dec         -6.072e-02  3.664e-02  -1.657   0.1285  
ow.Dec           3.553e-06  1.749e-06   2.032   0.0696 .
ow.frac.current  1.187e-01  4.807e-01   0.247   0.8100  
ice.Sep         -1.023e-05  8.558e-06  -1.195   0.2596  
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.09043 on 10 degrees of freedom
Multiple R-squared:  0.8579,	Adjusted R-squared:  0.7443 
F-statistic: 7.549 on 8 and 10 DF,  p-value: 0.002224

Hmmm… this looks messy.  There’s only one predictor with a slope significantly different from 0.  A high amount of variance in the original data is accounted for, but it smells like overfitting to me.  Let’s QC this model a bit by 1) checking for autocorrelations, 2) using AIC and relative likelihood to further eliminate predictors, and 3) selecting the final model with ANOVA.

## Check for autocorrelations using variance inflation factors

> vif(null.model)
         ao.Apr      ao.current         pdo.Mar         pdo.Aug         pdo.Dec          ow.Dec ow.frac.current 
       1.610203        1.464678        1.958998        2.854688        2.791749        1.772611        1.744653 

Surprisingly, all the parameters have acceptable vif scores (vif < 5).  Let’s proceed with relative likelihood.

## Define a function to evaluate relative likelihood.

rl <- function(aicmin, aici){
  return(exp((aicmin-aici) / 2))

## Construct multiple possible models, by adding parameters by order of abs(t value) in summary(null.lm).

model.1 <- lm(response ~ pdo.Mar, data = env.params)
model.2 <- lm(response ~ pdo.Mar + ao.current, data = env.params)
model.3 <- lm(response ~ pdo.Mar + ao.current + ow.Dec, data = env.params)
model.4 <- lm(response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec, data = env.params)
model.5 <- lm(response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec + ice.Sep, data = env.params)
model.6 <- lm(response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec + ice.Sep + ao.Apr, data = env.params)
model.7 <- lm(response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec + ice.Sep + ao.Apr + pdo.Aug, data = env.params)

## Collect AIC scores for models.

model.null.aic <- AIC(null.model)
model.1.aic <- AIC(model.1)
model.2.aic <- AIC(model.2)
model.3.aic <- AIC(model.3)
model.4.aic <- AIC(model.4)
model.5.aic <- AIC(model.5)
model.6.aic <- AIC(model.6)
model.7.aic <- AIC(model.7)

## Identify the model with the lowest AIC score.

> which.min(c(model.1.aic, model.2.aic, model.3.aic, model.4.aic, model.5.aic, model.6.aic, model.7.aic, model.null.aic))
[1] 5

So model.5 has the lowest AIC.  We need to check the relative likelihood of other models minimizing information loss. Models with values < 0.05 do not have a significant likelihood of minimizing information loss and can be discarded.

> rl(model.5.aic, model.1.aic)
[1] 8.501094e-05
> rl(model.5.aic, model.1.aic)
[1] 8.501094e-05
> rl(model.5.aic, model.2.aic)
[1] 0.00143064
> rl(model.5.aic, model.3.aic)
[1] 0.002415304
> rl(model.5.aic, model.4.aic)
[1] 0.7536875
> rl(model.5.aic, model.6.aic)
[1] 0.4747277
> rl(model.5.aic, model.7.aic)
[1] 0.209965
> rl(model.5.aic, model.null.aic)
[1] 0.08183346

Excellent, we can discard quite a few possible mode here; model.1, model.2, and model.3.  Last we’ll use ANOVA and the chi-squared test to see if any of the models are significantly different from model.4, the model with the fewest parameters.

> anova(model.4, model.5, test = 'Chisq')
Analysis of Variance Table

Model 1: response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec
Model 2: response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec + ice.Sep
  Res.Df      RSS Df Sum of Sq Pr(>Chi)
1     14 0.098626                      
2     13 0.086169  1  0.012457   0.1704

> anova(model.4, model.6, test = 'Chisq')
Analysis of Variance Table

Model 1: response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec
Model 2: response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec + ice.Sep + 
  Res.Df      RSS Df Sum of Sq Pr(>Chi)
1     14 0.098626                      
2     12 0.083887  2   0.01474   0.3485

> anova(model.4, model.7, test = 'Chisq')
Analysis of Variance Table

Model 1: response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec
Model 2: response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec + ice.Sep + 
    ao.Apr + pdo.Aug
  Res.Df      RSS Df Sum of Sq Pr(>Chi)
1     14 0.098626                      
2     11 0.082276  3   0.01635   0.5347

> anova(model.4, null.model, test = 'Chisq')
Analysis of Variance Table

Model 1: response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec
Model 2: response ~ ao.Apr + ao.current + pdo.Mar + pdo.Aug + pdo.Dec + 
    ow.Dec + ow.frac.current + ice.Sep
  Res.Df      RSS Df Sum of Sq Pr(>Chi)
1     14 0.098626                      
2     10 0.081778  4  0.016849   0.7247

Okay, there’s a whole lot going on there, they important thing is that none of the models are significantly different from the model with the fewest parameters. There are probably some small gains in the performance of those models, but at an increased risk of over fitting.  So model.4 is the winner, and we deem pdo.Mar, ao.current, ow.Dec, and pdo.Dec to be the best predictors of the response variable.  For those of you who are curious, that’s the March index for the Pacific Decadal Oscillation, the index for the Antarctic Oscillation when the cruise was taking place, within-pack ice open water extent in December (the month prior to the cruise), and the Pacific Decadal Oscillation index in December.  Let’s take a look at the model:

> summary(model.4)

lm(formula = response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec, 
    data = env.params)

      Min        1Q    Median        3Q       Max 
-0.177243 -0.037756 -0.003996  0.059606  0.114155 

              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.586e-02  5.893e-02  -0.609 0.552563    
pdo.Mar     -9.526e-02  2.208e-02  -4.315 0.000713 ***
ao.current  -3.127e-02  1.590e-02  -1.967 0.069308 .  
ow.Dec       4.516e-06  1.452e-06   3.111 0.007663 ** 
pdo.Dec     -8.264e-02  2.172e-02  -3.804 0.001936 ** 
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.08393 on 14 degrees of freedom
Multiple R-squared:  0.8287,	Adjusted R-squared:  0.7797 
F-statistic: 16.93 on 4 and 14 DF,  p-value: 2.947e-05

You’ll recall that the null model accounted for 74 % of the variance in the response variable, our final model accounts for 78 % and we’ve dropped several parameters.  Not bad!  Note that we never could have gotten to this point however, without the holistic search provided by glmnet.

I can still hear some grumbling about over fitting however, so let’s try to address that by bootstrapping.  We are crazy data-limited with only 19 observations, but let’s iteratively hold three observations back at random, build a model with the remaining 16 (using our identified parameters), and see how well each model predicts the three that were held back.  I created a function to do this so that I could repeat this exercise for lots of different response variables.

## The function requires as input the response vector, a vector of the predictors,
## the data frame of predictors, and the number of iterations you want to run.

bootstrap.model <- function(response.vector, predictors, data, n){
  predictor.df <-[,predictors])
  colnames(predictor.df) <- predictors
  model.predictions <- matrix(ncol = 3, nrow = 0)
  colnames(model.predictions) <- c('year', 'prediction', 'real')
  ## How many observations you need to withhold is determined by how many
  ## iterations you want to run.  For example, for 1000 iterations of
  ## unique random subsets of 19 observations, you need to withold 3.
  x <- 0
  boot <- 0
  for(y in (length(response.vector) - 1):1){
    while(x < n){
      boot <- boot + 1
      x <- y ** boot
  ## Now build n new models.
  for(i in 1:n){
    predict.i <-, boot)
    train.i <- which(!1:length(response.vector) %in% predict.i)
    temp.model <- lm(response.vector ~ ., data = predictor.df, subset = train.i) <-[predict.i,])
    colnames( <- predictors
    temp.predict <- predict(temp.model, newdata =, type = 'response')
    temp.actual <- response.vector[predict.i]
    temp.out <- matrix(ncol = 3, nrow = boot)
    try({temp.out[,1] <- names(response.vector)[predict.i]}, silent = T)
    temp.out[,2] <- temp.predict
    temp.out[,3] <- temp.actual
    model.predictions <- rbind(model.predictions, temp.out)
  mean.predicted <-[,2]), as.character(model.predictions[,3]), mean))
  sd.predicted <-[,2]), as.character(model.predictions[,3]), sd))
  ## Make an awesome plot.

  plot(mean.predicted[,1] ~ as.numeric(row.names(mean.predicted)),
     ylab = 'Predicted response',
     xlab = 'Observed response',
     pch = 19)
  for(r in row.names(mean.predicted)){
    lines(c(as.numeric(r), as.numeric(r)), c(mean.predicted[r,1], mean.predicted[r,1] + sd.predicted[r,1]))
    lines(c(as.numeric(r), as.numeric(r)), c(mean.predicted[r,1], mean.predicted[r,1] - sd.predicted[r,1]))
  abline(0, 1, lty = 2)
  mean.lm <- lm(mean.predicted[,1] ~ as.numeric(row.names(mean.predicted)))
  ## Capture stuff that might be useful later in a list.
  list.out <- list()
  list.out$model.predictions <- model.predictions
  list.out$mean.lm <- mean.lm
  list.out$mean.predicted <- mean.predicted
  list.out$sd.predicted <- sd.predicted

## And execute the function...

demo.boostrap <- bootstrap.model(response, c('pdo.Mar', 'ao.current', 'ow.Dec', 'pdo.Dec'), env.params, 1000)

Hey, that looks pretty good!  The dotted line is 1:1, while the solid line gives the slope of the regression between predicted and observed values.  The lines through each point give the standard deviation of the predictions for that point across all iterations.  There’s scatter here, but the model predicts low values for low observations, and high for high, so it’s a start…

Posted in Uncategorized | 2 Comments

Microbial community segmentation with R

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

## 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),
        = 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, = 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)

     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)

     main = '',
     type = "property",
     property = som.cluster.edges$cluster, = 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:

[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 | Leave a comment

New paper published in ISME Journal

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