Green Lab Notebook

So I've concatenated both MiSeq runs (not joining) and then clipping all at 250bp. The UPARSE documentation demands this, but it is not very straightforward in their pipeline. So I used FastX toolkit to do this.

So this script is a mix of UPARSE, QIIME, and FastX. This has resulted in keeping more sequences (especially in low yield samples) than any other method I tried.

Here is the R script regarding quality: http://rpubs.com/jfmeadow/27882

Mostly it came down to really low reverse read quality that tanked the total output number, as seen here:

Here is the workflow I am currently using to process sequence data using a combination of QIIME and UPARSE. Ann W put this together based on Mike Robeson's post here.

The only step that takes lots of time is the OTU table python script. No idea why it is so slow, and it might be worth rewriting that script to speed things up.

Definitions: 

We defined "counts" as the number of viable colonies growing on the plates after incubated for approximately 24 hours. We defined survival factor as the count divided by the mean count of the dark boxes for each organism. 

Results:

For E. coli, we found that there was a slight treatment effect for the UV boxes (viable colonies increased as exposure to light decreased). There was not a clear treatment effect in the visible boxes (viable colonies were relatively consistent throughout the box scheme). Note that the mean of the dark box counts was 11.73 colonies (we want between 20-200 colonies to try and decrease variation). The data for E. coli is shown below.

For Pseudomonas, we saw the viability of colonies drastically increase as UV light exposure decreased. We saw hardly any growth for the visible boxes, and the growth we did see was exposed to the least amount of visible light. The mean of the dark box counts was 37.55 colonies, which was in the target range of countable colonies (20-200). The data for Pseudomonas is shown below.

Consequences:

In this experiment, we saw that E. coli and Pseudomonas responded differently to the various light treatments. This was the first time that we tested E. coli on this scale. The data for Pseudomonas are consistent with what we saw in the most recent past trial. 

Next steps:

Github repository here: https://github.com/jfmeadow/ReproducibleDemo

Slides here: https://dl.dropboxusercontent.com/u/62653716/ReproDemo.pdf

This is a topic I've also been asked to talk about during my upcoming visit to San Luis Potosi, Mexico. So I'll get to reuse this potentially several times.

I tested my previous method by additionally comparing the slope of lines that are apparently parallel (fig below).

This is a very similar situation to the previous example, just a different data set. When I ran the old method on new data I got a very significant difference with a t-test comparing the means of these slopes, which I shouldn't be seeing. This poked a big hole in the way I had run things previously and as I thought about it, it made sense. The t-test is very sensitive to the n of the two populations being compared. So the 10,000 replicates from each bootstrapped sample was going to reveal a significant difference between the two almost no matter how similar the slopes were.

So I searched the internet a little to see what other folks had come up with. I found this stackexchange that provided an answer.

In a nutshell, what I ended up doing was to use the same bootstrapped slopes from my previous method, take the difference between the two populations of slopes (producing 10,000 differences), and then among those differences I calculated the proportion of those that were less than or equal to zero (which is basically the probability that the slopes were equal or overlapping). I then multiplied the proportion by two (to account for two tails) and then used the result to assess significance.

Here's the R code and output from what I did on the parallel lines above:

#take the differences, calculate the % of samples where the difference is <= 0 (where values are equal or overlap) and multiply that by 2. bootdiff.eea <- boot.eea.dry$t - boot.eea.wet$t hist(bootdiff.eea) overlap.eea <- sum(ifelse(bootdiff.eea <= 0, 1, 0)) # > overlap.eea # [1] 3352 2*overlap.eea/length(bootdiff.eea) # [1] 0.6704

Here's the histogram of differences:

So you can see that these two parallel slopes are definitely not significantly different.

Here's the same analysis with the data from the previous post:

bootdiff.otu <- boot.wet$t - boot.dry$t hist(bootdiff.otu) overlap.otu <- sum(ifelse(bootdiff.otu <= 0, 1, 0)) # > overlap.otu # [1] 1 2*overlap.otu/length(bootdiff.otu) # [1] 0.00020002

Histogram of differences:

There's a strong significant difference between these slopes, which is visually obvious.

Results:

All plates in the Visible boxes were clean. No viable colonies grew. UV plates displayed a nice shift of viability front-to-back based on expected light exposure levels. Here's the data plotted below. Black points are means, survival factor = (treatment count)/(mean of dark box counts).

Counts from the Dark boxes were fairly consistent within a box, but there were definitely differences between them. See this plot of the count data (black points are means, colors represent points from different boxes):

It looks like we're finally to the point where we've got apparently repeatable data between replicate boxes.

Next move:

I was working on the manuscript from my Master's work this weekend and had a statistical revelation regarding bootstrapping and comparing slopes that other lab folks thought would be worth sharing here.

A little context: the data I'm working with are from microbial community analysis of soil samples. I was working with a Bray-Curtis distance matrix and a geographic distance matrix to look at spatial patterns, via distance-decay.

So I wanted to perform distance-decay analysis on two subsets of my data: samples from wet sites and samples from dry sites. I generated a figure that basically looks like this:

Then I wanted to compare the slopes of these two lines and determine if the turnover of communities from dry sites was significantly different from wet sites. In order to do that I found a method in Horner-Devine's Nature paper from 2004 that utilizes a bootstrapping method. I do something slightly different, but I think achieve the same results.

A popular method to compare slopes is to use an ANCOVA, however because the points in my regressions aren't independent (they're pair-wise comparisons to generate a distance matrix) I can't do this. So instead, I bootstrapped the slopes of each line and then compared the bootstrapped slopes with a two sample t-test.

Here's my non-statistician way of explaining bootstrapping (feel free to correct me): you want to compare two statistics (like slopes of lines), but you can't just run a t-test on those two numbers because there's no sample variation to take into account which is necessary for statistical analysis. So you create a population of possible slopes by randomly subsampling from each data set -- this is bootstrapping. Then you get a population distribution that looks something like this:

Note that the slope from my original regression was 0.026.

I did this for both the dry and wet (or whatever categories you're comparing) slopes and then did a t-test to compare those two populations of slopes.

Here's snippets of the relevant code:

#define the statistic you want, in this case I want to bootstrap the slope of a regression slope <- function(formula, data, indices) {   d <- data[indices,] # allows boot to select sample fit <- lm(formula, data=d) return(summary(fit)$coefficients[2]) } #run the bootstrap boot.dry <- boot(data=table.dry, statistic=slope, R=9999, formula=log(X2)~log(X1)) boot.wet <- boot(data=table.wet, statistic=slope, R=9999, formula=log(X2)~log(X1)) #look at spread of values plot(boot.dry) #test to see if there's a difference between the two slopes t.test(boot.dry$t, boot.wet$t) #boom

These sites were really helpful for me:

Quick-R // Nonparametric bootstrapping: http://www.statmethods.net/advstats/bootstrapping.html

R in Ecology // Comparing two regression slopes by means of an ANCOVA: http://r-eco-evo.blogspot.com/2011/08/comparing-two-regression-slopes-by.html

Appendix to An R and S-PLUS Companion to Applied Regression // Bootstrapping Regression Models: http://cran.r-project.org/doc/contrib/Fox-Companion/appendix-bootstrapping.pdf

Here's a snippet of the code that ended up working:

<style type="text/css"> ul.mylist {list-style-position: inside;} ul.mylist li {text-indent: -1em; position: relative; left: 2em; margin-right: 2em;} ul.mylist li p {display: inline;} </style> <ul class="mylist"> <li>...

All samples were placed down the middle of the boxes. Location 1 is nearest the window, 5 is the furthest. Counts are diluted 10-7. Survival Factor is the proportion of viable counts from the light treatment sample compared to the average of the Dark box samples.

Raw E. coli counts:

Location Vis UV Dark 1 142 144 190 2 174 177 157 3 165 191 167 4 181 183 187 5 151 224 212 Average 162.6 183.8 182.6

E. coli Survival Factor:

Location Vis UV 1 0.78 0.79 2 0.95 0.97 3 0.90 1.05 4 0.99 1.00 5 0.83 1.23

Raw Pseudomonas counts:

Location Vis UV Dark 1 47 79 125 2 82 97 100 3 88 110 102 4 80 117 96 5 84 119 117 Average 76.2 104.4 108

Pseudomonas Survival Factor:

Location Vis UV 1 0.44 0.73 2 0.76 0.90 3 0.81 1.02 4 0.74 1.08 5 0.78 1.10

It's been a fantastic year so far.  Highlights

The script for the first chapter of our graphic novel about the urban microbiome, Noli timere, is complete.  Steve Green laid out the first storyboard - 24 pages, and is now working on the drawings.  All members of our team met at the Angouleme International Comics Festival, where I received my first "Auteur" Badge

Our J^4 Seagreass Microbiome project (Jonathan Eisen, Jenna Lang, Jay Stachowicz, Jess Green) is micrometers away from being formalized.  More on that soon.

Helene Morlon's lab has recently moved from Ecole Polytechnique to the Ecole Normal Superieure.  I will continue coming to Ecole, and also join Helene at ENS once a week.

I have given 7 of my 10 Blaise Pascal International Research Chair lectures, all of which I will be posting on this notebook in the coming weeks.

Installation instructions are here: http://joey711.github.io/phyloseq/install

So here is my code for

installing phyloseq

importing data

putting those data into normal R format for analysis.

Install:

They are not on CRAN, but they do have a nice bioconductor function for installation:

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

Upon installation, they will ask you to update lots of dependencies, and they have also made this totally painless. Just let them do it.

Import:

I downloaded the Jeon household bacteria data to play, and unzipped into a desktop folder.

directory <- '~/Desktop/study_2025_split_library_seqs_and_mapping/study_2025' biomfile <- paste0(directory, '_closed_reference_otu_table.biom') jeon.biom <- import_biom(biomfile, parseFunction=parse_taxonomy_greengenes) sdfile <- paste0(directory, '_mapping_file.txt') jeon.mapping <- import_qiime_sample_data(sdfile)

Now everything is in phyloseq format, which I find annoying, so I put into normal formatting.

Unpack:

jeon.table <- t(otu_table(jeon.biom)) jeon.map <- data.frame(jeon.mapping) jeon.taxa <- data.frame(tax_table(jeon.biom))

Then they are in easy formats, but need to be jived with one another.

jeon.map <- jeon.map[sort(row.names(jeon.map)), ] jeon.table <- jeon.table[sort(row.names(jeon.table)), ] identical(row.names(jeon.map), row.names(jeon.table)) #True identical(row.names(jeon.taxa), colnames(jeon.table)) #True

As those more detailed instructions say, you just need to extract and rerun the commands for the failed job in the file ending with _jobs.txt found in the output directory for the parallel script. There is also a script to help you figure out which job failed: identify_missing_files.py.

For example, this is part of a PBS script to rerun the failed job:

##execute program here: pick_otus.py -i picked_otus/prefilter_otus/POTU_B57u_/POTU_B57u_.10.fasta -r /home4/adamea/greengenes_13_5/gg_13_5_otus/rep_set/97_otus.fasta -m uclust_ref --suppress_new_clusters -o picked_otus/prefilter_otus/POTU_B57u_ -s 0.6 --max_accepts 20 --max_rejects 500 --stepwords 20 --w 12 ; mv picked_otus/prefilter_otus/POTU_B57u_/POTU_B57u_.10_otus.log picked_otus/prefilter_otus//POTU_B57u_.10_otus.log; mv picked_otus/prefilter_otus/POTU_B57u_/POTU_B57u_.10_otus.txt picked_otus/prefilter_otus//POTU_B57u_.10_otus.txt; mv picked_otus/prefilter_otus/POTU_B57u_/POTU_B57u_.10_failures.txt picked_otus/prefilter_otus//POTU_B57u_.10_failures.txt; mv picked_otus/prefilter_otus/POTU_B57u_/POTU_B57u_10_clusters.uc picked_otus/prefilter_otus//POTU_B57u_10_clusters.uc ; exit

Bioinformatics Pipeline

Raw Data -> OTU table & IDE

Using: ACISS, QIIME, FastX Toolkit

I received 137.9 GB of sequence data from the Genomics Core at the University of Oregon after submitting 36 samples for paired end 150 bp reads on the HiSeq run in MiSeq mode.

It was run with a companion of genomic DNA, due to the fact that a HiSeq can’t cope without sufficient sequence complexity. At first the sequence data was provided to me already demultiplexed, and this initially presented a problem because of the pipeline that I had originally vetted through QIIME.

I talked to Nick Stiffler, and he concluded the following: “The qiime docs show this when using files that were already demultiplexed, we need to create a special fasta label that includes the fields from the mapping file: http://qiime.org/documentation/file_formats.html-already-demultiplexed-samples. It doesn’t look like there is a script to handle this, but we could write one fairly easily.”

In the meantime, he gave me the raw data to use. In the future, I would have to use the QIIME script referenced above, along with a custom script to create a special fasta label.

It is important to remember that when logging on to the htseq server to download your data, you need to do so from your ACISS account, otherwise you will get an error: permission denied.

starting over with raw data (not demultiplexed)

Log into htseq server where the sequence data is, and download all relevant files.

scp [email protected]:/raid6p1/hiseq/ 140117_SN747_0364_BH8D14ADXX/ashley/Project_H8D14ADXX/ Sample_lane1/*.gz

Concatenate Read 1, Read 2 (barcode reads), and Read 3

cat *_R1_*.fastq < R1combined.fastq (Read 1) cat *_R2_*.fastq < R2combined.fastq (Barcodes) cat *_R3_*.fastq < R3combined.fastq (Read 3)

paste barcodes onto Read 1 sequence file

paste -d '\0' RawData_1.29.14/R2combined.fastq RawData_1.29.14/R1combined.fastq > PSG_R2R1.fastq

gets rid of double plusses

perl -i.bak -p -e 's/^\+\+\n/\+\n/' PSG_R2R1.fastq

We published the Lillis Dust paper in PLOS ONE this week, so I've also been fielding press for that. It has shown up in Gizmodo, Popular Science, Quartz, Fast Company, Futurity, and lots of other news sites that just printed the press release. So that's fun.

The surface paper was also accepted in Microbiome Journal so that will finally be published soon.