博文
宏基因组教程Metagenomics Tutorial (HUMAnN2)
||
之前我们介绍了microbiome_helper提供众多宏基因组学相关的实验、分析教程。
今天我们先详细讲解基于HUMAnN2的有参宏基因组分析流程 Metagenomics Tutorial (HUMAnN2) https://github.com/LangilleLab/microbiome_helper/wiki/Metagenomics-Tutorial-(Humann2))
此流程在我们之前学习的相关软件基础上，进一步完善了分析流程，包括数据质控、去宿主、
中间的物种和功能组成量我们之前有过介绍，
同时还提供了各步骤间的格式转换，以及下游的统计绘图。
说明：如果使用官方的虚拟机，可完全不用安装软件和修改原版代码即可运行测试数据。本文是基于我个人系统安装了最新版本软件，因此软件名称、安装位置、数据库位置都有变化，请用户根据自身情况选择学习。
分析流程
下载测试数据
mkdir -p ~/test/mh_meta17/v2 cd ~/test/mh_meta17/v2 wget http://kronos.pharmacology.dal.ca/public_files/tutorial_datasets/mgs_tutorial_Oct2017.zip unzip mgs_tutorial_Oct2017.zip cd mgs_tutorial_Oct2017 gunzip raw_data/*
了解输入文件
# 实验设计 less map.txt # 样品数 wc -l map.txt # 测序数据 head -n 8 raw_data/p136C_R1.fastq
输入文件1：实验设计
主要包括样品名，分组(如正常/癌症)，相关属性
Sample Id Sample Type Sex Individual p136C Cancer M p136 p136N Normal M p136 p143C Cancer M p143
输入文件2：测序文件，每个样品1个或1对fastq/fq文件
@SRR3586062.883556 CTTGGGGCTGCTGAGCTTCATGCTCCCCTCCTGCCTCAAGGACAATAAGGAGATCTTCGACAAGCCTGCAGCAGCTCGCATCGACGCCCTCATCGCTGAGG + CCCFFFFFHHHHHIJJJJJJJIJIJJJJGIJDGIJEIIJIJJJJJJJJIJJJJIJJIJJJJJHHHFFFFECEEEDDDDD?BDDDDDDBDDDDDDDDBBBDD
软件安装和环境变量
export PATH=$PATH:/conda2/bin
详细安装方法见文末软件安装部分。为什么我不把安装目录永久添加在环境变量，因为conda也会影响环境变量，导致我其它工具依赖关系出错，因此特殊流程可以临时添加更安全，自己清楚在每个流程在哪里，即使添加也要放在$PATH后面才安全。
不可能存在一个环境满足所有的依赖关系，所以近几年conda, docker会非常流行。推荐安装首选conda，如果仍搞不定，再用docker(小提示，每个bioconda软件都提供docker下载，只是使用更麻烦但成功率更高)。
序列质控和去宿主
# 新版本脚本名称变更，指定bowtie2数据库、trimmomatic位置 kneaddata -h # 显示帮助 kneaddata -i raw_data_example/p144C_R1.fastq -i raw_data_example/p144C_R2.fastq -o kneaddata_out -v \ -db ~/ref/host/Homo_sapiens \ --trimmomatic /mnt/bai/public/bin/Trimmomatic-0.36/ --trimmomatic-options "SLIDINGWINDOW:4:20 MINLEN:50" \ -t 9 --bowtie2-options "--very-sensitive --dovetail" --remove-intermediate-output
-h 显示帮助
-v 显示运行信息，方便观察运行进展和出错找原因
-i 为输入fastq文件，双端需输两次
-o 输出结果目录
-db 指定bowtie2索引
—trimmomatic 指定质控程序位置，默认为/usr/bin/trimmomatic-0.36.jar
—trimmomatic-options 质控选项，4碱基滑窗，质量大于20，最小长度50
-t 设置线程数，不要超过9
—bowtie2-options 比对选项
—remove-intermediate-output 清理中间文件
批处理质控
你不可能实验只一个文件，一般最小2个实验组，每组三次重复，需要批量处理
这里使用parallel命令，依赖map文件中样品名，或文件列表对测序结果并行处理。注意不同版本下参数可能不同，如处理多参数时，原文虚拟机20170322版本为—links参数，而我的系统中20141022版本为—xapply参数
# 批量读取成对文件作为输入
parallel -j 3 --xapply 'echo {1} {2}' ::: raw_data/*_R1.fastq ::: raw_data/*_R2.fastq
# 批处理样品
parallel -j 3 --xapply 'kneaddata -i {1} -i {2} -o kneaddata_out -v \
-db ~/ref/host/Homo_sapiens \
--trimmomatic /mnt/bai/public/bin/Trimmomatic-0.36/ --trimmomatic-options "SLIDINGWINDOW:4:20 MINLEN:50" \
-t 9 --bowtie2-options "--very-sensitive --dovetail" --remove-intermediate-output' \
::: raw_data/*_R1.fastq ::: raw_data/*_R2.fastq质控后结果统计
kneaddata_read_count_table --input kneaddata_out --output kneaddata_read_counts.txt
合并双端
# 清理宿主污染至指定目录备用
mkdir kneaddata_out/contam_seq
mv kneaddata_out/*_contam*fastq kneaddata_out/contam_seq
# 本质上只合并了1/2
concat_paired_end.pl -p 9 --no_R_match -o cat_reads kneaddata_out/*_paired_*.fastq
# 可选方法2：我感觉是合并4个文件，应该重命名序列ID，不则双端序列重名字，shell命令合并单样品
for i in `tail -n+2 map.txt | cut -f 1`;do \
cat kneaddata_out/${i}_R1_kneaddata_paired_1.fastq kneaddata_out/${i}_R1_kneaddata_paired_2.fastq kneaddata_out/${i}_R1_kneaddata_unmatched_1.fastq kneaddata_out/${i}_R1_kneaddata_unmatched_2.fastq | awk '{if(NR%4==1) print "@"NR; else print $0}' > cat_reads/${i}.fastq; done计算功能和代谢通路
humman2是一套分析流程，它包括调用metaphlan流程来分析物种组成，和自身分析功能基因和代谢通路组成。安装方法见本文软件安装部分。
humann2 --threads 9 --input cat_reads/p144C.fastq --output humann2_out/
多样品批量并行处理，-j参数请尽量小于你你电脑核心数(线程数的一半)，否则效率反而会低且系统卡顿。
parallel -j 9 'humann2 --threads 9 --input {} --output humann2_out/{/.}' ::: cat_reads/*fastq多样品物种和功能组成合并为矩阵/表
merge_metaphlan_tables.py precalculated/metaphlan2_out/*tsv > metaphlan2_merged.tsv # 转换为spf格式方便stamp分析 metaphlan_to_stamp.pl metaphlan2_merged.tsv > metaphlan2_merged.spf # 删除样品名中多余字符，和株水平行 sed -i 's/_metaphlan_bugs_list//g' metaphlan2_merged.spf # 删除株行：新数据库新增行 sed -i '/t__[^\t]*\t/d' metaphlan2_merged.spf
STAMP软件统计绘图
可以对上述结果使用STAMP打开，鼠标点点可完成常见统计分析的可视化。主页 http://kiwi.cs.dal.ca/Software/STAMP
可参考我们之前的教程
可进一步在软件主页学习第一手英文教程。
整理humann2功能结果
humann2_join_tables --input precalculated/humann2_out/ --file_name pathabundance --output humann2_pathabundance.tsv # 标准化为百分比 humann2_renorm_table --input humann2_pathabundance.tsv --units relab --output humann2_pathabundance_relab.tsv # 分层结果 humann2_split_stratified_table --input humann2_pathabundance_relab.tsv --output ./ # 筛选某个通路结果 head -n 1 humann2_pathabundance_relab_unstratified.tsv >> humann2_pathabundance_relab_LACTOSECAT-PWY.tsv grep "LACTOSECAT-PWY" humann2_pathabundance_relab_unstratified.tsv >> humann2_pathabundance_relab_LACTOSECAT-PWY.tsv # 格式化为stamp输入 sed 's/_Abundance//g' humann2_pathabundance_relab_LACTOSECAT-PWY.tsv > humann2_pathabundance_relab_LACTOSECAT-PWY.spf sed -i 's/# Pathway/MetaCyc_pathway/' humann2_pathabundance_relab_LACTOSECAT-PWY.spf
使用STAMP分析结果
软件安装
本流程虽然只有简单的几步，但是背后却是几十步，需要的软件极多，也许安装和配置的复杂程序，反倒是分析流程难度的几倍。
作者提供了Virtalbox虚拟机 https://github.com/LangilleLab/microbiome_helper/wiki/Microbiome-Helper-Virtual-Box，有20GB大小，可以直接学习，但灵活性不高，运行效率也低。有服务器的朋友应该更喜欢正常安装软件使用更高效方便。
依赖软件清单如下，这里仅对我使用的系统没有的软件进行安装说明，其它软件请读者自行查看软件官方帮助安装。
https://github.com/LangilleLab/microbiome_helper/wiki/Requirements
流程相关脚本
为完成流程各软件间的衔接，需要写很多脚本，作者整理了几十个常用脚本工具，可下载并添加环境
cd ~/github git clone git@github.com:LangilleLab/microbiome_helper.git echo 'export PATH=$PATH:~/github/microbiome_helper' >> ~/.bashrc # 安装perl模块 cpan install File::Basename Getopt::Long List::Util Parallel::ForkManager Pod::Usage
质控和去宿主kneadData
https://bitbucket.org/biobakery/kneaddata/wiki/Home
KneadData是一款宏基因组和宏转录组测序数据质控的流程，其主要功能包括使用Trimmomatic序列质控，bowtie2比对至宿主和PhiX基因组去除宿主和测序平衡序列。
#以下三个方式均不可用 #pip install kneaddata # 无权限 #pip install kneaddata --user # 不满足依赖关系 #pip3 install kneaddata --user # 安装成功，位于~/.local/lib/python3.5/site-packages/kneaddata，配置了依赖软件和数据库，在bowtie步仍报错 # /conda2安装 /conda2/bin/conda install kneaddata # 列出主要命令 ll /conda2/bin/kneaddata* echo 'export PATH=$PATH:/conda2/bin' >> ~/.bashrc # 或者临时添加conda2加入环境变量 export PATH=$PATH:/conda2/bin # 如果你还不可用，使用docker吧，https://quay.io/repository/biocontainers/kneaddata，每个conda都有docker可下载，但也需要更大权限
查看可用数据
kneaddata_database
下载人类和小鼠基因组数据库
mkdir -p ~/ref/host cd ~/ref/host kneaddata_database --download human_genome bowtie2 ./ kneaddata_database --download mouse_C57BL bowtie2 ./
创建定制数据库
本质上就是把指定的基因组建一个bowtie2的索引
https://bitbucket.org/biobakery/kneaddata/wiki/Home#markdown-header-create-a-custom-database
以下载的EMBL plants上的拟南芥基因组为例
# bowtie2-build <reference> <db-name> cd ~/ref/ath bowtie2-build Arabidopsis_thaliana.TAIR10.dna.toplevel.fa bt2ebi
humann2安装
详见之前的文章
比对安装在/conda/bin中，如下命令可以Ubuntu 16.04中添加环境变量
echo 'export PATH=$PATH:/conda/bin' >> ~/.bashrc
猜你喜欢
10000+：肠道细菌 人体上的生命 宝宝与猫狗 梅毒狂想曲 提DNA发Nature 实验分析谁对结果影响大 Cell微生物专刊
写在后面
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外1500+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。
学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组”
点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA
https://wap.sciencenet.cn/blog-3334560-1115288.html
上一篇：微生物组助手——最易学的扩增子、宏基因组分析流程
下一篇：STAMP：扩增子、宏基因组统计分析神器(中文帮助文档)
