分析流程
- 下载测试数据
- 了解输入文件
- 软件安装和环境变量
- 序列质控和去宿主
- 质控后结果统计
- 合并双端
- 计算功能和代谢通路
- 多样品物种和功能组成合并为矩阵/表
- STAMP软件统计绘图
- 整理humann2功能结果
软件安装
- 流程相关脚本
- 质控和去宿主kneadData
- humann2安装
- 猜你喜欢
- 写在后面

之前我们介绍了microbiome_helper提供众多宏基因组学相关的实验、分析教程。

微生物组助手——最易学的扩增子、宏基因组分析流程

今天我们先详细讲解基于HUMAnN2的有参宏基因组分析流程 Metagenomics Tutorial (HUMAnN2) https://github.com/LangilleLab/microbiome_helper/wiki/Metagenomics-Tutorial-(Humann2)

此流程在我们之前学习的相关软件基础上，进一步完善了分析流程，包括数据质控、去宿主、

中间的物种和功能组成量我们之前有过介绍，
- MetaPhlAn2一条命令获得宏基因组物种组成
- HUMAnN2：人类微生物组统一代谢网络分析2

同时还提供了各步骤间的格式转换，以及下游的统计绘图。

说明：如果使用官方的虚拟机，可完全不用安装软件和修改原版代码即可运行测试数据。本文是基于我个人系统安装了最新版本软件，因此软件名称、安装位置、数据库位置都有变化，请用户根据自身情况选择学习。

分析流程

下载测试数据

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

可参考我们之前的教程
- 微生物组间差异分析神器-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安装

详见之前的文章
- HUMAnN2：人类微生物组统一代谢网络分析2

比对安装在/conda/bin中，如下命令可以Ubuntu 16.04中添加环境变量

echo 'export PATH=$PATH:/conda/bin' >> ~/.bashrc

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外1500+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA

宏基因组教程Metagenomics Tutorial (HUMAnN2)相关推荐

宏基因组有参分析和无参分析差异
宏基因组有参分析和无参分析差异分析流程解决问题结果差异宏基因组流程综述本文参考宏基因组教程Metagenomics Tutorial (HUMAnN2) 分析流程有参流程:质控–物种组成 ...
Nature子刊：HUMAnN2实现宏基因组和宏转录组种水平功能组成分析
HUMAnN2实现宏基因组和宏转录组种水平功能组成分析 Species-level functional profiling of metagenomes and metatranscriptomes ...
Nature Method：HUMAnN2实现宏基因组和宏转录组种水平功能组成分析
文章目录 HUMAnN2实现宏基因组和宏转录组种水平功能组成分析简介热心肠日报导读摘要主要图表图1. HUMAnN2分层式搜索在同类软件中准确率最高图2. 人类核心微生物组的贡献多样性图 ...
你想要的宏基因组-微生物组知识全在这(190101)
文章目录征稿.转载.合作文章分类导航目录精选文章推荐培训.会议.招聘科研经验软件和数据库扩增子分析宏基因组分析 Linux与Shell R统计绘图实验设计与技术基础知识一作解读 ...
你想要的宏基因组-微生物组知识全在这(181101)
文章目录征稿.转载.合作文章分类导航目录精选文章推荐培训.会议.招聘广告科研经验软件和数据库扩增子分析宏基因组分析 R统计绘图实验设计与技术基础知识必读综述高分文章套路解读科 ...
你想要的宏基因组-微生物组知识全在这(181001)
文章目录征稿.转载.合作文章分类导航目录精选文章推荐培训.会议.招聘广告科研经验软件和数据库扩增子分析宏基因组分析 R统计绘图实验设计与技术基础知识必读综述高分文章套路解读科 ...
你想要的宏基因组-微生物组知识全在这(180801)
本文"宏基因组"原创,更多文章点我跳转公众号阅读宏基因组/微生物组是当今世界科研最热门的研究领域之一,为加强本领域的技术交流与传播,推动中国微生物组计划发展,中科院青年科研人员创 ...
你想要的宏基因组-微生物组知识全在这(180601)
征稿.转载.合作文章分类导航目录精选文章推荐培训.会议.招聘广告科研经验软件和数据库宏基因组分析 R统计绘图实验设计与技术基础知识必读综述高分文章套路解读科普视频-寓教于乐写在 ...
你想要的宏基因组-微生物组知识全在这(180901)
征稿.转载.合作文章分类导航目录精选文章推荐培训.会议.招聘广告科研经验软件和数据库扩增子分析宏基因组分析 R统计绘图实验设计与技术基础知识必读综述高分文章套路解读科普视频-寓 ...

宏基因组教程Metagenomics Tutorial (HUMAnN2)