查看原文
其他

宏基因组教程Metagenomics Tutorial (HUMAnN2)

metagenome 宏基因组 2019-07-15

之前我们介绍了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安装

详见之前的文章

猜你喜欢

10000+:肠道细菌 人体上的生命 宝宝与猫狗 梅毒狂想曲 提DNA发Nature 实验分析谁对结果影响大  Cell微生物专刊

系列教程:微生物组入门 Biostar 微生物组  宏基因组

专业技能:生信宝典 学术图表 高分文章 不可或缺的人

一文读懂:宏基因组 寄生虫益处 进化树

必备技能:提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析:图表解读 分析流程 统计绘图

16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun

在线工具:16S预测培养基 生信绘图

科研经验:云笔记  云协作 公众号

编程模板 Shell  R Perl

生物科普  生命大跃进  细胞暗战 人体奥秘  

写在后面

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

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

点击阅读原文,跳转最新文章目录阅读

    您可能也对以下帖子感兴趣

    文章有问题?点此查看未经处理的缓存