CSHL Computational Genomics - November, 2011 -- Metagenomics I

  1. This workshop will run from the Unix/MacOS command line.
    Students on Mac's can login to their workstations, and start the Terminal application (in Applications/Utilitles.
    Students on PC's need to login to courses.cshl.edu.

  2. Once in a terminal window, run the script /ecg/seqprg/scripts/init_meta1.sh. Type:
    If things have worked properly, you should have the directory meta1_work in your home directory, and it should contain several files.
    cd meta1_work

  3. To do the analysis, we will run the mothur program to analyze microbial communities. We are using an example from http://www.mothur.org/wiki/Schloss_SOP.

    To run the mothur program, type:

    All commands in mothur look like functions, and need to end with (), for example: help() or quit().

  4. In this workshop, we will be doing microbial community analysis on a sample set of 16S rRNA sequences from a human stool sample. The mothur program will allow us to address the following questions:
    1. What is the taxonomic makeup of the sample.
    2. How diverse is the community (what is the dynamic range of abundance); species "evenness", species "richness".
    3. How to compare two microbial community samples.

  5. In real analyses, the sequences must be pre-processed to remove bar codes, primers, and non-rRNA sequences. This pre-processing has already been done for the sequences labeled "final.*". The preprocessing steps are listed here

  6. The first step classifies the sequences taxonomically, and bins them into clades. This has already been done, producing the files final.names,final.taxonomy,final.group,final.fasta, etc. These are linked into your meta1_work directory. In addition, a time-consuming step produces final.dist.

    To characterize the taxonomic makeup, we first cluster sequences into OTUs (taxa, clades), use make.shared to count sequence abundance within those OTUs.

    The lines below show commands that can be copied and pasted into the command line of mothur. Lines beginning with # are comments to explain the steps; you only need to copy lines that do not begin with #.

    #Cluster sequences into OTU -- we are using the command cluster.split to do that as it allows us to cluster
    # sequences according a taxonomic level like order or family
    cluster.split(fasta=final.fasta, taxonomy=final.taxonomy, name=final.names, taxlevel=3, processors=4)
    # The make.shared creates a file that represent the number of times
    # that an OTU is observed in multiple samples
    make.shared(list=final.an.list, group=final.groups, label=0.03)
    # Since some samples might have better coverage, sub-sample to get a dataset
    # with the same number of sequences per sample
    # assign a consensus taxonomy to each OTU
    classify.otu(list=final.an.list, name=final.names, taxonomy=final.taxonomy, label=0.03, cutoff=80)
    # Calculate the relative abundance of taxa
    phylotype(taxonomy=final.taxonomy, name=final.names, label=1)
    # The make.shared creates a file that represent the number of times that an OTU is observed in multiple samples
    make.shared(list=final.tx.list, group=final.groups, label=1)
    # Since some samples might have better coverage, subsample to get a dataset
    # with the same number of sequences per sample
    sub.sample(shared=final.tx.shared, size=400)
    # assign a consensus taxonomy to each OTU
    classify.otu(list=final.tx.list, name=final.names, taxonomy=final.taxonomy, label=1)
  7. Once the sequences have been classified, we can build a tree that represents the population:
    #Here we want to build a phylogenetic tree
    sub.sample(fasta=final.fasta, name=final.names, group=final.groups, persample=T, size=400)
    dist.seqs(fasta=final.subsample.unique.fasta, output=lt, processors=2)
    #Build Tree
    You can visualize the tree with TreeView (in /ecg/Applications/TreeView).

  8. To quantify population diversity, collector's curves are produced:
    #Generates collector's curves
    collect.single(shared=final.an.0.03.subsample.shared, calc=chao-invsimpson, freq=100)

  9. To compare two populations, one can use a heat map or venn diagram:
    #Generate heatmaps and a venn diagram to compare samples
    heatmap.bin(scale=log2, numotu=50) 
    #PCOA Analysis to Compare samples
    dist.shared(shared=final.an.0.03.subsample.shared, calc=thetayc-jclass)
    #Non-metric multidimensional scaling 

    These programs create .svg files, which you should be able to visualize with a web browser.

Course Home Page