PICRUSt2 Streamlined Tutorial (v2.1.3 beta) - picrust/picrust2 GitHub Wiki

This is the streamlined version of the PICRUSt2 tutorial, which focuses on running commands. You can find the full tutorial here. You can see descriptions of the arguments and options used by each script by typing the script name followed by --help. All past and present official PICRUSt2 tutorials are listed here. Make sure you are looking at the correct version of the tutorial.

PICRUSt2 is only available for Mac OS X and Linux systems (not for Windows unless you are running a Linux emulator). Before running the installation commands below you will need to install Miniconda.

To run this tutorial you will need at least 16 GB of RAM, however you can still run through most of this tutorial with less RAM. The key memory-intensive step is the read placement, which you can avoid by simply downloading the expected output file and proceeding with the tutorial. When running your data you can run the QIIME2 plugin of PICRUSt2 and use SEPP for read placement. Using this approach is much less memory intensive, but currently there are fewer options available for the QIIME2 plugin.

The easiest way to install PICRUSt2 is with bioconda, which will install a new environment which contains all the software needed to run PICRUSt2.

conda create -n picrust2 -c bioconda -c conda-forge picrust2=2.1.3_b

When using conda you need to specifically activate environments to use them:

conda activate picrust2

First, download, unzip, and enter the folder of tutorial files:

wget http://kronos.pharmacology.dal.ca/public_files/picrust/picrust2_tutorial_files/chemerin_16S.zip
unzip chemerin_16S.zip
cd chemerin_16S

You should see the three files below when you type ls -1:

metadata.tsv
seqs.fna
table.biom

These files are processed datafiles rather than raw reads. We can take a look at the input FASTA file with less (note that with less you can go back to the terminal by typing q).

less seqs.fna

You should see:

>86681c8e9c64b6683071cf185f9e3419
CGGGCTCAAACGGGGAGTGACTCAGGCAGAGACGCCTGTTTCCCACGGGACACTCCCCGAGGTGCTGCATGGTTGTCGTC
AGCTCGTGCCGTGAGGTGTCGGCTTAAGTGCCATAACGAGCGCAACCCCCGCGTGCAGTTGCTAACAGATAACGCTGAGG
ACTCTGCACGGACTGCCGGCGCAAGCCGCGAGGAAGGCGGGGATGACGTCAAATCAGCACGGCCCTTACGTCCGGGGCGA
CACACGTGTTACAATGGGTGGTACAGCGGGAAGCCAGGCGGCGACGCCGAGCGGAACCCGAAATCCACTCTCAGTTCGGA
TCGGAGTCTGCAACCCGACTCCGTGAAGCTGGATTCGCTAGTAATCGCGCATCAGCCATGGCGCGGTGAATACGTTCCCG
>b7b02ffee9c9862f822c7fa87f055090
AGGTCTTGACATCCAGTGCAAACCTAAGAGATTAGGTGTTCCCTTCGGGGACGCTGAGACAGGTGGTGCATGGCTGTCGT
CAGCTCGTGTCGTGAGATGTTGGGTTAAGTCCCGCAACGAGCGCAACCCTTGTCATTAGTTGCCATCATTAAGTTGGGCA
CTCTAATGAGACTGCCGGTGACAAACCGGAGGAAGGTGGGGATGACGTCAAGTCATCATGCCCCTTATGACCTGGGCTAC
ACACGTGCTACAATGGACGGTACAACGAGAAGCGAACCTGCGAAGGCAAGCGGATCTCTTAAAGCCGTTCTCAGTTCGGA
CTGTAGGCTGCAACTCGCCTACACGAAGCTGGAATCGCTAGTAATCGCGGATCAGCACGCCGCGGTGAATACGTTCCCGG
...

This file contains the Amplicon Sequence Variants (ASVs) of 16S rRNA gene sequences found across the mouse stool samples. The id of each ASV is written on the lines starting with ">", so the id of the first sequence is "86681c8e9c64b6683071cf185f9e3419".

The next file is a bit tricky to view since it is a BIOM file which is binary encoded (not easily viewed with a plain-text viewer). We can get a sense of what the first few columns and rows look like by using the following command:

biom head -i table.biom

You should see:

# Constructed from biom file
#OTU ID	100CHE6KO	101CHE6WT	102CHE6WT	103CHE6KO	104CHE6KO
86681c8e9c64b6683071cf185f9e3419	716.0	699.0	408.0	580.0	958.0
b7b02ffee9c9862f822c7fa87f055090	51.0	98.0	96.0	91.0	223.0
663baa27cddeac43b08f428e96650bf5	112.0	29.0	45.0	35.0	75.0
abd40160e03332710641ee74b1b42816	0.0	0.0	0.0	0.0	0.0
bac85f511e578b92b4a16264719bdb2d	56.0	4.0	63.0	131.0	90.0

Notice, that the first column ("#OTU ID") contains the ASV ids in the FASTA file above. The first column of this table needs to match the ids in the FASTA file. The additional columns represent different samples with the counts representing the number of reads within each of those samples. Importantly, the input table should contain read counts rather than relative abundance. You can input a rarefied table, but that is not required.

We can get more information about the rest of the table using this command:

biom summarize-table -i table.biom

The easiest way to run PICRUSt2 is with the picrust2_pipeline.py script, which you can see here. This script automatically will run all of the steps that are described below. However, for the purposes of the tutorial we will run each step individually so that it is easier to follow.

First make and enter a new folder for writing the PICRUSt2 output files.

mkdir picrust2_out_pipeline
cd picrust2_out_pipeline

When running the below commands you can swap in the number of available cores for the 1 in -p 1, which will allow many steps to be run in parallel and speed up the pipeline.

The first step of PICRUSt2 is to insert your study ASVs into a reference tree (details).

You can run this step with this command (1 minute and 13.8 GB of RAM on 1 processor):

place_seqs.py -s ../seqs.fna -o out.tre -p 1 \
              --intermediate intermediate/place_seqs

The required input is the FASTA of the ASV sequences. The key output file is out.tre, which is a tree in newick format of your study ASVs and reference 16S sequences. If you do not have enough RAM to run this command locally you can simply download the expected out.tre output file here.

Run hsp.py, which will take 4 min and 1 GB of RAM on 1 processor with this dataset. Note that if you see a "FutureWarning" you can ignore that message.

hsp.py -i 16S -t out.tre -o marker_predicted_and_nsti.tsv -p 1 -n

hsp.py -i EC -t out.tre -o EC_predicted.tsv -p 1

The output files of these commands are marker_predicted_and_nsti.tsv (download) and EC_predicted.tsv (download).

The metagenome_pipeline.py script will output tables of gene family predictions per sample (see details), which should run in ~3 seconds:

metagenome_pipeline.py -i ../table.biom -m marker_predicted_and_nsti.tsv -f EC_predicted.tsv \
                       -o EC_metagenome_out --strat_out 

In the above command the --strat_out option indicates that we also want the stratified output file (the ASVs contributing to each function). There are several output files created within the EC_metagenome_out directory:

• EC_metagenome_out/pred_metagenome_unstrat.tsv (download) - overall EC number abundances per sample.
• EC_metagenome_out/pred_metagenome_strat.tsv (download) - EC number abundances per sample stratified by the contributing ASVs.
• EC_metagenome_out/seqtab_norm.tsv (download) - the ASV abundance table normalized by predicted 16S copy number.
• EC_metagenome_out/weighted_nsti.tsv (download) - the mean NSTI value per sample (when taking into account the relative abundance of the ASVs). This file can be useful for identifying outlier samples in your dataset. In PICRUSt1 weighted NSTI values < 0.06 and > 0.15 were suggested as good and high, respectively. The cut-offs can be useful for getting a ball-park of how your samples compare to other datasets, but a weighted NSTI score > 0.15 does not necessarily mean that the predictions are meaningless.

The last major step of the PICRUSt2 pipeline is to infer pathway-level abundances with pathway_pipeline.py (details), which you can run with this command (< 2 min):

pathway_pipeline.py -i EC_metagenome_out/pred_metagenome_strat.tsv \
                    -o pathways_out -p 1

The output of this script are in the pathways_out folder. These include both unstratified (download) and stratified (download) MetaCyc pathway abundances, similar to the EC number tables above. Note that this stratified output was output since a stratified gene family table was input.

Finally, it can be useful to have a description of each functional id in the output abundance tables. The below commands will add these descriptions as new column in gene family and pathway abundance tables (details):

add_descriptions.py -i EC_metagenome_out/pred_metagenome_unstrat.tsv -m EC \
                    -o EC_metagenome_out/pred_metagenome_unstrat_descrip.tsv

add_descriptions.py -i pathways_out/path_abun_unstrat.tsv -m METACYC \
                    -o pathways_out/path_abun_unstrat_descrip.tsv

Note that the mapping of function ids to descriptions can also be found in the files within this folder of the PICRUSt2 repository: picrust2/picrust2/default_files/description_mapfiles/.

As mentioned above, there are many possible ways to analyze the PICRUSt2 output. STAMP is one tool that can be used that requires no background in scripting languages. This tool is very useful for visualizing microbiome data, but note that there are better alternatives for conducting statistical tests with compositional data (such as the ALDEx2 R package).

1: Install STAMP (this has to be installed on a local machine, not a server since it has a graphical interface)

• Windows: Download and install.
• Mac or Linux:

conda deactivate
conda install -c bioconda stamp

Then simply type the following to start:

stamp

2: If you ran PICRUSt2 on a remote server you will need to download path_abun_unstrat_descrip.tsv and metadata.tsv to your own computer. You can use WinSCP (for Windows) or CyberDuck (for Mac).

3: Start up STAMP and load in the path_abun_unstrat_descrip.tsv as the Profile File, and the metadata.tsv file as the Group metadata file. You can now explore the data overall with heatmaps and ordination tools like PCA or alternatively compare boxplots of the relative abundance of individual predicted pathways.