PICRUST2 16s扩增子群落功能预测 安装和使用

测试系统：WSL2 Ubuntu 20.04LTS
1.照着Github一键安装

conda create -n picrust2 -c bioconda -c conda-forge picrust2=2.4.1 -y
conda activate picrust2

错误可能如下：
（1）卡在Solving environment，安装不上去，多半网络问题，多执行几次就好
尽管我使用的也是清华源，但是因为前期设置的.condarc没有使用新版conda的参数设定，所以默认使用的频道依旧为conda官方（海外地址）的，因此solve起来特别缓慢，之后直接使用conda create -n picrust2 picrust2=2.4.1 -y解决
我原先的清华源设置：

show_channel_urls: true
channels:
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge/
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
  - defaults

现在清华源推荐的设置：

show_channel_urls: true
channels:
  - defaults
default_channels:
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main
  - https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
custom_channels:
  conda-forge: https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud
  bioconda: https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud

（2）安装完成，执行picrust2_pipeline.py -s study_seqs.fna -i study_seqs.biom -o picrust2_out_pipeline -p 1总是在metagenome_pipeline.py这一步报错，未知问题

2.照着Github分步安装

wget https://github.com/picrust/picrust2/archive/v2.4.1.tar.gz
tar xvzf v2.4.1.tar.gz
cd picrust2-2.4.1/

conda env create -f  picrust2-env.yaml
conda activate picrust2
pip install --editable .

pytest

测试也总是出错，不能通过place_seqs.py这一脚本

3.针对1中（2）和2解决手段：
更换系统，使用Ubuntu 18.04LTS一键安装，执行下来不会报错，记录该方法在此，目前尚不清楚是何种原因导致Ubuntu 20.04LTS执行失败

picrust2_pipeline.py -s study_seqs.fna -i study_seqs.biom -o picrust2_out_pipeline -p 16

之后添加注释

cd picrust2_out_pipeline
add_descriptions.py -i EC_metagenome_out/pred_metagenome_unstrat.tsv.gz -m EC -o EC_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz

add_descriptions.py -i KO_metagenome_out/pred_metagenome_unstrat.tsv.gz -m KO  -o KO_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz

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

最后对KEGG通路进行分级（脚本参考自YongxinLiu/EasyMicrobiome）

os.makedirs("./KEGG_classification")
cd KEGG_classification
wget -c https://raw.githubusercontent.com/YongxinLiu/EasyMicrobiome/main/script/summarizeAbundance.py
wget -c https://raw.githubusercontent.com/YongxinLiu/EasyMicrobiome/main/kegg/KO1-4.txt
zcat ../KO_metagenome_out/pred_metagenome_unstrat.tsv.gz > KEGG.KO.txt
python3 summarizeAbundance.py -i KEGG.KO.txt -m KO1-4.txt -c 2,3,4 -s ',+,+,' -n raw -o KEGG

实际执行步骤分解：

①将扩增子序列变体（或OTU）放入参考系统发育

PICRUSt2包装HMMER以将研究序列放入参考多序列比对中，然后使用EPA-NG或SEPP将这些序列放入参考系统发育中。
所谓的“研究序列”将是典型工作流程下的代表性OTU和/或ASV。
工具GAPPA用于将生成的.jplace对象转换为newick格式。
Tips：您的输入研究序列需要在正链上！

place_seqs.py -s study_seqs.fna -o placed_seqs.tre -p 16 \
              --intermediate placement_working

②运行隐藏状态预测，得到每个代表序列预测的16S拷贝数、E.C.数（EnzymeCommission）和KO丰度（KEGG Orthology）。

PICRUSt2包装了castor R包以运行隐藏状态预测(hsp)来预测基因家族丰度。

hsp.py -i 16S -t placed_seqs.tre -o marker_nsti_predicted.tsv.gz -p 16 -n
#请注意，NSTI 值将添加到 16S 预测表中（因为指定了参数-n）
hsp.py -i EC -t placed_seqs.tre -o EC_predicted.tsv.gz -p 16
hsp.py -i KO -t placed_seqs.tre -o KO_predicted.tsv.gz -p 16

③预测测序样本中EC和KO基因家族丰度（通过16S序列丰度调整基因家族丰度）

脚本metagenome_pipe.py读取序列丰度表(以bio、TSV或mother shared文件格式的OTUs或asv的丰度)、预测标记基因丰度文件和预测基因家族丰度文件(最后两个文件由hsp.py输出)。
序列丰度应该是读数计数，而不是相对丰度。
它将根据预测的标记基因数量对输入的序列丰度表进行归一化。然后，它将确定每个样本的预测功能状况。
Tips：即使输入文件是BIOM格式，输出文件也以制表符分隔。

metagenome_pipeline.py -i study_seqs.biom \
                       -m marker_nsti_predicted.tsv.gz \
                       -f EC_predicted.tsv.gz \
                       -o EC_metagenome_out

metagenome_pipeline.py -i study_seqs.biom \
                       -m marker_nsti_predicted.tsv.gz \
                       -f KO_predicted.tsv.gz \
                       -o KO_metagenome_out

④根据预测的EC数量丰度推断MetaCyc通路丰度和覆盖率

通路丰度使用与HUMAnN2相同的方法计算，基于与途径内反应相关的基因家族丰度(默认将E.C.数重新分组为MetaCyc反应途径)。
可以输入结构化或非结构化通路映射文件（默认情况下映射文件是结构化的），这将根据必需基因家族的存在来识别可能存在的通路集。
pathway_pipeline.py脚本使用两个默认映射文件。这些文件是默认指定的，因此您无需自己指定！但是，了解此脚本默认执行的操作很有用。
首先使用此映射文件default_files/pathway_mapfiles/ec_level4_to_metacyc_rxn.tsv将EC号重新分组为MetaCyc RXNs；
然后，这些MetaCyc RXNs可以使用这个映射文件default_files/pathway_mapfiles/metacyc_path2rxn_struc_filt_pro.txt来推断MetaCyc通路丰度。第二个映射文件包含对原核生物中发现的MetaCyc通路子集的通路的反应图。

pathway_pipeline.py -i EC_metagenome_out/pred_metagenome_unstrat.tsv.gz \
                    -o pathways_out \
                    --intermediate pathways_working \
                    -p 16

⑤在基因家族和通路丰度表中添加描述作为新列

add_descriptions.py是一个方便的脚本，它将在您的基因家族或通路丰度表中添加一列，对应于每个功能类别的快速描述。这些描述在picrust2/default_files/description_mapfiles。
您还可以使用自定义映射文件。

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

add_descriptions.py -i KO_metagenome_out/pred_metagenome_unstrat.tsv.gz -m KO -o KO_metagenome_out/pred_metagenome_unstrat_descrip.tsv.gz

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