#jamesmeadow

[greenlabnotebook-blog]

After lots of variations on this analysis pipeline, I've found that joining paired ends unnecessarily removes lots of sequences. Likely this is due to quality dropoff at about 250bp (fig below).

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.

############# commands to run interactive node on aciss cluster qsub -q fatnodes -I module load usearch module load fastx_toolkit/0.0.13 module load qiime/1.8.0 printenv PATH # check to make sure lots of qiime dependencies are loaded. # Qiime has so many dependencies, it is difficult to load them all. # So might have to reload module version. This makes no sense. source /usr/local/packages/Modules/setmodule 3.2.9 source /usr/local/packages/Modules/setmodule 3.2.10 ################ ##### Check quality - this goes into R script Pickle2014/R/quality/ # split_libraries_fastq.py -v -q 0 -i raw1/r1readTRIM.fastq -b raw1/barcodesRenamed.fastq -o splitLib1F/ -m map.txt --barcode_type 16 # split_libraries_fastq.py -v -q 0 -i raw2/r1readTRIM.fastq -b raw2/barcodesRenamed.fastq -o splitLib2F/ -m map.txt --barcode_type 16 # split_libraries_fastq.py -v -q 0 -i raw1/r2readTRIM.fastq -b raw1/barcodesRenamed.fastq -o splitLib1R/ -m map.txt --barcode_type 16 # split_libraries_fastq.py -v -q 0 -i raw2/r2readTRIM.fastq -b raw2/barcodesRenamed.fastq -o splitLib2R/ -m map.txt --barcode_type 16 # Ended up combining forward reads and just analyzing these instead of joining. cat raw1/r1readTRIM.fastq raw2/r1readTRIM.fastq > seqs.fastq cat raw1/barcodesRenamed.fastq raw2/barcodesRenamed.fastq > barcodes.fastq split_libraries_fastq.py -v -q 0 --store_demultiplexed_fastq -i seqs.fastq -b barcodes.fastq -o splitLib/ -m map.txt --barcode_type 16 # -n 300 # trim to 250 length. This is not straightforward or well documented in UPARSE, so # farm out to fastx. fastx_trimmer -l 250 -i splitLib/seqs.fastq -o splitLib/seqs.trimmed.fastq -Q33 # get quality stats usearch -fastq_stats splitLib/seqs.trimmed.fastq -log splitLib/seqs.stats.log # remove low quality reads mkdir qF usearch -fastq_filter splitLib/seqs.trimmed.fastq -fastq_maxee 0.5 -fastaout qF/seqs.filtered.fasta # dereplicate sequences. Last step with files separate. mkdir deRep usearch -derep_fulllength qF/seqs.filtered.fasta -output deRep/seqs.filtered.derep.fasta -sizeout # filter singletons - This rids sigletons - Decided to do without # mkdir filterSingles # usearch -sortbysize deRep/seqs.filtered.derep.fasta -minsize 2 -output filterSingles/seqs.filtered.derep.mc2.fasta # clusterOTUs mkdir OTUs usearch -cluster_otus deRep/seqs.filtered.derep.fasta -otus OTUs/seqs.filtered.derep.repset.fasta # reference chimera check mkdir chiCheck usearch -uchime_ref OTUs/seqs.filtered.derep.repset.fasta -db scripts/gold.fa -strand plus -nonchimeras chiCheck/seqs.filtered.derep.repset.nochimeras.fasta # label OTUs using puthon script from UPARSE mkdir labelOTUs python scripts/fasta_number.py chiCheck/seqs.filtered.derep.repset.nochimeras.fasta OTU_ > labelOTUs/seqs.filtered.derep.repset.nochimeras.otus.fasta # match original quality filtered reads back to otus - this is with bash derep workaround. mkdir matchOTUs usearch -usearch_global qF/seqs.filtered.fasta -db labelOTUs/seqs.filtered.derep.repset.nochimeras.otus.fasta -strand plus -id 0.97 -uc matchOTUs/otu.map.uc # make otu table mkdir otuTable python scripts/uc2otutab_mod.py matchOTUs/otu.map.uc > otu-table.txt # convert to biom biom convert --table-type="OTU table" -i otu-table.txt -o otu-table.biom # **use QIIME 1.7, not 1.8** Dependency problem # assign taxonomy assign_taxonomy.py -t gg_13_5_otus/taxonomy/97_otu_taxonomy.txt -r gg_13_5_otus/rep_set/97_otus.fasta -i labelOTUs/seqs.filtered.derep.repset.nochimeras.otus.fasta -o assigned_taxonomy # add taxonomy to BIOM table biom add-metadata --sc-separated taxonomy --observation-header OTUID,taxonomy --observation-metadata-fp assigned_taxonomy/seqs.filtered.derep.repset.nochimeras.otus_tax_assignments.txt -i otu-table.biom -o otu_table.biom # check sequencing depth. # print_biom_table_summary.py -i otu_table.biom ## for qiime <1.8 biom summarize-table -i otu_table.biom -o otu_table_summary.txt # for qiime >=1.8

[greenlabnotebook-blog]

After reviewing the joined read results for pickle box analysis, I decided to join the forward reads from both runs, and avoid joining altogether. The reverse reads were pretty poor across the board, so we were losing way to many sequences and joining others incorrectly.

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:

The x and y axes here are equal at the top sequence count. Thus the forward read in the second run had the most sequences, and by combining both forward runs, I can get the most from the total number of sequences.

[greenlabnotebook-blog]

I've been trying to find alternative methods to the standard QIIME OTU clustering. UPARSE (published 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.

# Split libraries with QIIME split_libraries_fastq.py -v -q 0 --store_demultiplexed_fastq -i $COMBINED/seqs.fastq -b $COMBINED/barcodes.fastq -o splitLib/ -m map.txt --barcode_type 16 # -n 300 # get quality stats usearch -fastq_stats splitLib/seqs.fastq -log splitLib/seqs.stats.log # remove low quality reads - trimmed short seqs - presumeably didn"t join correctly. mkdir qF usearch -fastq_filter splitLib/seqs.fastq -fastq_maxee 0.5 -fastaout qF/seqs.filtered.tmp.fasta -fastq_minlen 400 # -fastq_trunclen 296 sed 's/>/>barcodelabel=/' qF/seqs.filtered.tmp.fasta > qF/seqs.filtered.fasta # dereplicate sequences. Last step with files separate. mkdir deRep usearch -derep_fulllength qF/seqs.filtered.fasta -output deRep/seqs.filtered.derep.fasta -sizeout # filter singletons mkdir filterSingles usearch -sortbysize deRep/seqs.filtered.derep.fasta -minsize 2 -output filterSingles/seqs.filtered.derep.mc2.fasta # clusterOTUs mkdir OTUs usearch -cluster_otus filterSingles/seqs.filtered.derep.mc2.fasta -otus OTUs/seqs.filtered.derep.mc2.repset.fasta # reference chimera check mkdir chiCheck usearch -uchime_ref OTUs/seqs.filtered.derep.mc2.repset.fasta -db scripts/gold.fa -strand plus -nonchimeras chiCheck/seqs.filtered.derep.mc2.repset.nochimeras.fasta # label OTUs using puthon script from UPARSE mkdir labelOTUs python scripts/fasta_number.py chiCheck/seqs.filtered.derep.mc2.repset.nochimeras.fasta OTU_ > labelOTUs/seqs.filtered.derep.mc2.repset.nochimeras.otus.fasta # match original quality filtered reads back to otus - this is with bash derep workaround. mkdir matchOTUs usearch -usearch_global qF/seqs.filtered.fasta -db labelOTUs/seqs.filtered.derep.mc2.repset.nochimeras.otus.fasta -strand plus -id 0.97 -uc matchOTUs/otu.map.uc # make otu table mkdir otuTable # python scripts/uc2otutab.py matchOTUs/otu.map.uc > otuTable/seqs.filtered.derep.mc2.repset.nochimeras.otu-table.txt python scripts/uc2otutab_jl.py matchOTUs/otu.map.uc > otuTable/seqs.filtered.derep.mc2.repset.nochimeras.otu-table.txt #### still slow - running pbs script Monday afternoon. # convert to biom biom convert --table-type="OTU table" -i otuTable/seqs.filtered.derep.mc2.repset.nochimeras.otu-table.txt -o otuTable/seqs.filtered.derep.mc2.repset.nochimeras.otu-table.biom # assign taxonomy assign_taxonomy.py -t gg_13_5_otus/taxonomy/97_otu_taxonomy.txt -r gg_13_8_otus/rep_set/97_otus.fasta -i labelOTUs/seqs.filtered.derep.mc2.repset.nochimeras.otus.fasta -o assigned_taxonomy # add taxonomy to BIOM table biom add-metadata --sc-separated taxonomy --observation-header OTUID,taxonomy --observation-metadata-fp assigned_taxonomy/seqs.filtered.derep.mc2.repset.nochimeras.OTUs_tax_assignments.txt -i otuTable/seqs.filtered.derep.mc2.repset.nochimeras.otu-table.biom -o otuTable/otu_table.biom

[greenlabnotebook-blog]

We received great feedback on the knitr documents that were submitted with the Lillis microbial surface paper that was published in Microbiome Journal. The editors actually wrote a great piece about our efforts! So I've been asked to speak recently about how to pull this off. This has given me an opportunity to create some teaching materials. Mostly composed of a small subset of those data, along with analysis scripts and example manuscript documents all created in the dynamic analysis document style, using knitr and pandoc.

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.

Lately we've been working really hard getting the next round of PickleBox studies off the ground, including IRB and planning with ESBL. We're also almost ready to submit our phone microbiome paper.

[greenlabnotebook-blog]

Finally dug into the phyloseq package in R. This is designed to make microbiome studies easier to analyze in R. I found that most of the functionality is sort of useless for me - since I already have functions and routines for everything they do. But they do have a nice quick function to read OTU tables in .biom format, and a similarly easy qiime mapping file entry.

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

So now everything is in normal ecology format for R analysis!

[greenlabnotebook-blog]

I've been working all week to migrate lots of data onto our new storage system. We're trying to maintain a consistent directory structure for each project so that we can all navigate easily in each others' projects well after we each leave.

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.

I also submitted the pickle box paper last week.

[greenlabnotebook-blog]

Finished a draft of the picklebox manuscript today and handed it off for editing! It is coming along really nicely and the study gets better and better as I write and dig into results and figures. Rare.

As for manuscript writing in markdown, it is going really well. Except tables are a real bummer. Converting tables between markdown, latex, and word is a hassle, and editing tables in anything other than excel or R is also a hassle. So for now, I'm keeping tables out (in their own excel file) until the final version. Then they can be plugged in pretty easily. Other than that, manuscript writing in markdown is wonderful and pretty much hassle free!

I'll be absent from the lab for a couple of weeks. Enjoy!

[greenlabnotebook-blog]

Converting figures for journal submission is always a beast. Often it requires clunky proprietary software, and so I've been looking for an alternative. I found a really easy solution on Linux, and it will also work on a Mac.

On Ubuntu Linux, ImageMagick comes pre-installed, and contains the convert function for lots of different formats. It is also really easy to install on Mac if port is already installed (port is a command-line install tool that makes installation more Linux-like). [EDIT] To instal port on a mac, you also have to have xcode installed, which is somewhat of a hassle, but worth it in the long run.

To install ImageMagick on mac:

sudo port install ImageMagick

Then, to convert a figure:

convert <to> <from>

For instance, most of my figures are vector (.pdf) but need to be submitted as .png since it is a lossless raster format:

convert image.pdf image.png

This is low resolution, so change dpi (density per inch)

convert -density 500 image.pdf image.png

This also works for tiff, postscript, jpeg, etc.

[EDIT] Adam found this list of convert arguments: http://www.imagemagick.org/script/convert.php