CBW 2016 PICRUSt tutorial

==Introduction==

This lab will walk you through the basic steps of using PICRUSt to make functional predictions (e.g. predicted metagenome) from your 16S data.

It uses an OTU table that has already been generated for use with PICRUSt. See detailed instructions on how to do this using closed reference picking [http://picrust.github.io/picrust/tutorials/otu_picking.html here] or using open-reference picking [http://github.com/mlangill/microbiome_helper#picrust-workflow-for-16s-data here]

The data we will be using in this lab comes from the stool of three groups of mice that are of different ages (e.g. young, middle, and old).

== Preliminaries ==

=== Work directory === Create a new directory that will store all of the files created in this lab.

rm -rf ~/workspace/module_picrust
mkdir -p ~/workspace/module_picrust
cd ~/workspace/module_picrust

Now we need to download the data using 'wget':

wget https://www.dropbox.com/sh/a35f90j8eh3r23j/AADzQ9zLrEud5xHAHG8kKxlua?dl=1 -O picrust_data.zip

Now decompress the data using "unzip" command:

unzip picrust_data.zip
rm picrust_data.zip

== Main Lab Steps ==

The first step is to correct the OTU table based on the predicted 16S copy number for each organism in the OTU table:

normalize_by_copy_number.py -i otus.biom -o otus_corrected.biom

Note that this is just a normal OTU table which then could be used with all of the other tools you already learned today.

If you want to look at the before and after correction you can use the biom tools to convert it to plain text:

biom convert -i otus_corrected.biom -o otus_corrected.txt --to-tsv --header-key taxonomy
biom convert -i otus.biom -o otus.txt --to-tsv --header-key taxonomy

Now you can look at them using less:

less otus.txt

# Constructed from biom file
#OTU ID	9Y-June1	10Y-June1	8Y-May23	10Y-May23	6Y-June1	9Y-May23	Y7-Aug14	Y7-Aug15	6Y-May23	M11-Aug14	M11-Aug15	M11-Jul13	11M-May23	M13-Jul13	13M-May23	2E-Aug14	2E-Aug15	2E-May24	4E-June1	1E-Aug16	1E-May23	taxonomy
181348	1.0	0.0	0.0	0.0	0.0	1.0	0.0	0.0	1.0	0.0	0.0	0.0	0.0	0.0	0.0	6.0	15.0	3.0	4.0	7.0	0.0	k__Bacteria;  p__Firmicutes;  c__Clostridia;  o__Clostridiales;  f__Lachnospiraceae;  g__Coprococcus;  s__
318732	0.0	0.0	1.0	0.0	0.0	0.0	0.0	0.0	2.0	5.0	9.0	7.0	1.0	5.0	3.0	0.0	2.0	0.0	0.0	0.0	0.0	k__Bacteria;  p__Firmicutes;  c__Clostridia;  o__Clostridiales;  f__;  g__;  s__
244484	0.0	0.0	0.0	2.0	0.0	1.0	0.0	1.0	4.0	0.0	2.0	0.0	2.0	1.0	0.0	0.0	1.0	0.0	2.0	0.0	1.0	k__Bacteria;  p__Firmicutes;  c__Clostridia;  o__Clostridiales;  f__Ruminococcaceae;  g__;  s__
(etc.)

less otus_corrected.txt

#OTU ID	9Y-June1	10Y-June1	8Y-May23	10Y-May23	6Y-June1	9Y-May23	Y7-Aug14	Y7-Aug15	6Y-May23	M11-Aug14	M11-Aug15	M11-Jul13	11M-May23	M13-Jul13	13M-May23	2E-Aug14	2E-Aug15	2E-May24	4E-June1	1E-Aug16	1E-May23	taxonomy
181348	0.333333333333	0.0	0.0	0.0	0.0	0.333333333333	0.0	0.0	0.333333333333	0.0	0.0	0.0	0.0	0.0	0.0	2.0	5.0	1.0	1.33333333333	2.33333333333	0.0	k__Bacteria;  p__Firmicutes;  c__Clostridia;  o__Clostridiales;  f__Lachnospiraceae;  g__Coprococcus;  s__
318732	0.0	0.0	0.333333333333	0.0	0.0	0.0	0.0	0.0	0.666666666667	1.66666666667	3.0	2.33333333333	0.333333333333	1.66666666667	1.0	0.0	0.666666666667	0.0	0.0	0.0	0.0	k__Bacteria;  p__Firmicutes;  c__Clostridia;  o__Clostridiales;  f__;  g__;  s__
244484	0.0	0.0	0.0	1.0	0.0	0.5	0.0	0.5	2.0	0.0	1.0	0.0	1.0	0.5	0.0	0.0	0.5	0.0	1.0	0.0	0.5	k__Bacteria;  p__Firmicutes;  c__Clostridia;  o__Clostridiales;  f__Ruminococcaceae;  g__;  s__
(etc.)

Ok, no lets actually make our functional predictions of KEGG Ortholog (KOs) predictions using the corrected OTU table as input:

predict_metagenomes.py -i otus_corrected.biom -o ko_predictions.biom

We can check out these KO predictions again by converting the BIOM file first:

biom convert -i ko_predictions.biom -o ko_predictions.txt --to-tsv --header-key KEGG_Description    

# Constructed from biom file
#OTU ID	9Y-June1	10Y-June1	8Y-May23	10Y-May23	6Y-June1	9Y-May23	Y7-Aug14	Y7-Aug15	6Y-May23	M1
1-Aug14	M11-Aug15	M11-Jul13	11M-May23	M13-Jul13	13M-May23	2E-Aug14	2E-Aug15	2E-May24	4E-June1	1E
-Aug16	1E-May23	KEGG_Description
K01365	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.
0	0.0	0.0	cathepsin L [EC:3.4.22.15]
K01364	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.
0	0.0	0.0	streptopain [EC:3.4.22.10]
K01361	18.0	20.0	9.0	4.0	11.0	4.0	9.0	6.0	6.0	7.0	10.0	11.0	9.0	11.0	8.0	32.0	8.0	15.0	17
.0	8.0	9.0	lactocepin [EC:3.4.21.96]
K01360	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.
0	0.0	0.0	proprotein convertase subtilisin/kexin type 2 [EC:3.4.21.94]
K01362	3587.0	3559.0	3868.0	3428.0	3872.0	3462.0	3432.0	1913.0	2436.0	3219.0	3248.0	3081.0	3372.0	2602.0	3494.0	3566.0	3527.0	2616.0	31
33.0	3457.0	2212.0	None
K02249	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.
0	0.0	0.0	competence protein ComGG
K05841	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.
0	0.0	0.0	sterol 3beta-glucosyltransferase [EC:2.4.1.173]

Default predictions from PICRUSt are KOs (KEGG Orthologs) but PICRUSt can also predict COGs and Rfams.

PICRUSt can also collapse KOs to KEGG Pathways. Note that one KO can map to many KEGG Pathways so a simple mapping wouldn't work here. Instead, we use the PICRUSt script "categorize_by_function.py":

categorize_by_function.py -i ko_predictions.biom -c KEGG_Pathways -l 3 -o pathway_predictions.biom

Again lets look at the output:

biom convert -i pathway_predictions.biom -o pathway_predictions.txt --to-tsv --header-key KEGG_Pathway

# Constructed from biom file
#OTU ID	9Y-June1	10Y-June1	8Y-May23	10Y-May23	6Y-June1	9Y-May23	Y7-Aug14	Y7-Aug15	6Y-May23	M1
1-Aug14	M11-Aug15	M11-Jul13	11M-May23	M13-Jul13	13M-May23	2E-Aug14	2E-Aug15	2E-May24	4E-June1	1E
-Aug16	1E-May23	KEGG_Pathways
1,1,1-Trichloro-2,2-bis(4-chlorophenyl)ethane (DDT) degradation	11.0	21.0	10.0	7.0	14.0	4.0	8.0	1.0	4.0	0.0	0.0	0.
0	1.0	0.0	0.0	1.0	1.0	2.0	4.0	2.0	1.0	Metabolism; Xenobiotics Biodegradation and Metabolism; 1,1,1-Trichloro-2,2
-bis(4-chlorophenyl)ethane (DDT) degradation
ABC transporters	200982.0	174898.0	195247.0	255298.0	147766.0	254328.0	306731.0	490225.0	36
3852.0	217743.0	231867.0	239470.0	201328.0	237358.0	189880.0	199125.0	342119.0	294970.0	21
3939.0	229000.0	451627.0	Environmental Information Processing; Membrane Transport; ABC transporters
Adherens junction	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.
0	0.0	0.0	0.0	0.0	Cellular Processes; Cell Communication; Adherens junction
Adipocytokine signaling pathway	6486.0	6300.0	7408.0	6562.0	7205.0	6982.0	6139.0	4365.0	5299.0	7160.0	7528.0	6977.0	8475.0	6064.0	7827.0	74
04.0	7462.0	6411.0	7082.0	7654.0	5580.0	Organismal Systems; Endocrine System; Adipocytokine signaling pathway
African trypanosomiasis	28.0	25.0	40.0	42.0	23.0	26.0	188.0	33.0	22.0	62.0	63.0	43.0	19.0	12.0	29.0	19.0	24
.0	12.0	9.0	22.0	9.0	Human Diseases; Infectious Diseases; African trypanosomiasis
Alanine, aspartate and glutamate metabolism	94807.0	90632.0	103163.0	103640.0	96543.0	104717.0	106172.0	112557.0	10
5979.0	93152.0	100320.0	98573.0	101380.0	90366.0	98759.0	100108.0	113079.0	103468.0	98339.0	104441.0	115040.0	Metabolism; Amino Acid Metabolism; Alanine, aspartate and glutamate metabolism

PICRUSt can directly connect the OTUs that are contributing to each KO by using the ''metagenome_contributions.py'' script:

metagenome_contributions.py -i otus_corrected.biom -l K01727,K01194,K01216,K11049,K00389,K00449 -o metagenome_contributions.txt

Now lets look at getting more detail for the individual KOs that we focused on with the metagenome_contributions.py command from a few steps ago. You can browse the file using 'less':

less metagenome_contributions.tab

The output should look like this:

Gene Sample OTU GeneCountPerGenome OTUAbundanceInSample CountContributedByOTU ContributionPercentOfSample ContributionPercentOfAllSamples K01727 9Y-June1 190026 1.0 1.66666666667 1.66666666667 0.251889168766 0.000792700810933 K01727 9Y-June1 4331760 3.0 1.0 3.0 0.453400503778 0.00142686145968 K01727 9Y-June1 2594570 1.0 0.333333333333 0.333333333333 0.0503778337531 0.000158540162187 K01727 9Y-June1 1106050 1.0 0.333333333333 0.333333333333 0.0503778337531 0.000158540162187 K01727 9Y-June1 3090117 1.0 0.2 0.2 0.0302267002519 9.5124097312e-05 K01727 9Y-June1 1051299 1.0 0.75 0.75 0.113350125945 0.00035671536492 K01727 9Y-June1 2617854 1.0 0.333333333333 0.333333333333 0.0503778337531 0.000158540162187

Each line in this file relates how much a single OTU (third column) contributes to a single KO (first column) within a single sample (second column). The fifth column contains the actual relative abundance contributed by this OTU, and the other columns contain other information about the abundance of the OTU the percentage contribution of this OTU.

You could use your favourite plotting program (e.g. excel, sigmaplot, etc) to plot the information from columns 1-3 and column 5. As an example of what the output might look, I have created the following image:

File:K00449 genus.png

This plot shows that the large increase in K00449 within sample 25 is contributed by the genus ''Citrobacter''.

You can use STAMP with the corrected OTU table by first converting it using the Microbiome Helper script:

biom_to_stamp.py -m taxonomy otus_corrected.biom > otus_corrected.spf

Then we use a script from Microbiome Helper to convert the BIOM file into a STAMP profile file:

biom_to_stamp.py -m KEGG_Pathways pathway_predictions.biom > pathway_predictions.spf

Now download the pathway_predictions.spf file and the map.tsv file to your local computer and load these files within STAMP (File->Load).

Change Profile Level to "Level 3" which actually represents the KEGG Pathways. Then change the "Group Field" (top right) to "Age_Group".

Apply a multiple test correction and then view the individual KEGG Pathways using a "Box Plot" (instead of PCA plot). What is the most significant KEGG Pathway?

If you like you can explore other visualizations with STAMP or attempt to load KEGG Modules or KOs instead within STAMP.