宏基因组单个样本数据处理流程笔记

前言
数据预处理
- 质量控制
- 去除接头序列
- 去除宿主序列
物种注释
- Kraken2注释
- Krona制图
序列组装
- MEGAHIT序列组装
- SPAdes序列组装
功能注释
- HUMAnN功能注释
- RGI耐药注释
- VFDB毒力注释
后记

前言

衷心感谢这个开源互助的网络时代，能够让我在36岁高龄重新回到程序员的行列。虽然人生会因为各种原因起起落落落落落落，但还好家学渊源不至于无事可做。我在落落落落落落期间接触生物信息基础，竟然发现：这不就是当年奥赛做错的动态规划大题么？真是山不转水转……在此特别感谢陈同和刘永鑫两位老师的生物信息交流平台，以及GitHub上各个开源工具的作者和用户们的无私分享。
本篇博文将总结分享我在进行宏基因组数据处理中的一些流程和经验，并根据工作进度和涉猎范围进行完善补充，尤其欢迎各位朋友能够在留言区进行交流。本文不单独介绍各个软件的安装调试过程，而是集中在具体使用。硬件使用的是自费购买搭建的Intel NUC10i7FNH平台，系统是Ubuntu 20.04LTS版，环境控制是Anaconda。有Dell的运算服务器当然更好，但有也不等于就一定能用，没有也不等于不能开工。代码只要想写，就一定能写，这跟让不让你做手术是不太一样的。

数据预处理

跟若干家公司打过交道，要么拒绝提供数据，要么仅仅提供处理过的CleanData。为了能够确保流程和数据的可靠性，特别是对于以后要撰写论文过程中的严谨性，我要求提供下机时生成的RawData，以及Adaptor和Barcode等配套数据。
这一部分内容主要参考羊城迷鹿的文章宏基因组测序流程，涉及数据质量控制、去除接头和去除宿主几个步骤。

质量控制

使用的是FastQC进行质量控制，使用也很方便：

# 这里会直接在当前目录下生成html报告文件，--threads 12参数表示使用12个线程。
$ fastqc --threads 12 seqs.fastq.gz

其实对任何生成的fastq文件，无论是RawData、CleanData还是中间处理过程中生成的文件，都最好能够随时进行质控，做到心中有数。

去除接头序列

目前我还没有真实的执行过去除接头序列的步骤，因为前面几家公司都没有提供这些原始数据供验证。我目前备用的工具是trimmomatic，等我实际操作之后再进行补充。

去除宿主序列

我使用bwa去除宿主序列。之前用过别的一些软件，但感觉bwa去除得干净一些。比对用的数据库是人类基因组的hg38数据库，貌似这个应该是最新的版本了。
注意在使用之前需要先对下载的数据库文件hg38.fa建立索引。

$ bwa index -a bwtsw ~/hg38/hg38.fa

建立索引之后就可以开始去除宿主序列，基本思路就是先进行数据比对，并将比对上的序列进行反转去除，剩下未能比对成功的，就是去除宿主之后的目标序列：

# 使用12个线程将fastq文件在hg38数据库中进行比较，并转换成sam文件。
$ bwa mem -t 12 -M ~/hg38/hg38 seqs.fastq.gz > seqs.sam
# 将sam文件转换成bam，方便后续进行处理。
$ samtools view -bS seqs.sam > seqs.bam
# 筛选出没有在hg38当中被比中的序列，这里的具体原理可以看看软件的帮助文件或者羊城迷鹿的博文。
# 这里-f使用的4，羊城迷鹿的原文是12，是因为我用的是single序列，而对方用的是paired序列，所以根据软件的使用要求数值有差异。
$ samtools view -b -f 4 -F 256 seqs.bam > seqs.hostfree.bam
# 对bam文件内容进行排序。
$ samtools sort -n seqs.hostfree.bam -o seqs.hostfree.sorted.bam
# 对排序的结果转换为序列文件。
$ bedtools bamtofastq -i seqs.hostfree.sorted.bam -fq seqs.hostfree.fastq

获得的seqs.hostfree.fastq文件会用于后续各个步骤。需要注意的是，上述步骤未必能够保证宿主序列能够被完全去除，后续的物种注释当中仍然会出现少量的宿主序列，但比例已经少了很多。

物种注释

物种注释是宏基因组分析当中需要首先解决的问题：是什么，或者谁在那里。

Kraken2注释

物种注释的方法很多，但是我个人喜欢结果比较丰富的Kraken2。这个软件安装过程非常波折，因为要用非常巨大的物种数据库建立哈希表，下载这个数据库就花了很长时间。这个软件会在运行时把这个硕大的哈希表整个载入内存，所以内存小了无法运行。我的NUC最大安装内存是64G，所以还需要建立足够的虚拟内存。大家可以搜索一下具体的方法，主要是使用Linux的dd命令，随后用swapon载入。我建立了一个150G的虚拟空间，基本够用。
Kraken2的使用方法也是很简单：

# 这里使用我放在~/kraken_db/目录的数据库，12个线程全部工作。
$ kraken2 --db ~/kraken_db/ --threads 12 seqs.hostfree.fastq --output seqs.hostfree.kra --report seqs.hostfree.tsv

随着主机风扇呼啸，鼠标键盘卡到几乎凝固。一切结束之后会出现两个文件：–output选项的seqs.hostfree.kra是比对的报告文件，如果没有–output这个选项，这个报告会被直接扔到屏幕上，比弹幕护体更加疯狂；–report选项里的内容是我们需要的，包括比对出界门纲目科属种各个层级的物种丰度。这个格式看起来还是有点儿吃力的，其实还可以改成MetaPhlAn的输出格式。这些不同的输出格式适用于不同的后续处理流程。

# 在--report选项之后加入--use-mpa-style选项，就可以进行MetaPhlAn的格式输出。
$ kraken2 --db ~/kraken_db/ --threads 12 seqs.hostfree.fastq --output seqs.hostfree.kra --report seqs.hostfree.mpa.tsv --use-mpa-style

seqs.hostfree.tsv这种格式用于绘制后续的Krona花冠图。seqs.hostfree.mpa.tsv主要用于后续多个样本的数据分析。seqs.hostfree.kra这个比对报告文件并非无用，也可以用于花冠图的绘制。

Krona制图

绘制花冠图当然要用到Krona，但是在使用之前需要对前面步骤获得的文件进行处理，便于该软件识别绘图。这里用到的工具是Kraken的贡献者之一Jennifer Lu所制作的KrakenTools工具包，其中有一个工具是kreport2krona.py。把这个文件拷贝到工作目录，GitHub上有非常清楚的应用范例，具体使用方法如下：

# 特别注意，这里输入的文件是seqs.hostfree.tsv，不是seqs.hostfree.mpa.tsv。
$ python kreport2krona.py -r seqs.hostfree.tsv -o seqs.hostfree.kro

获得的seqs.hostfree.kro就可以拿来输入Krona生成花冠图：

# 这里务必注意命令的大小写。
$ ktImportText seqs.hostfree.kro -o seqs.hostfree.html

点开seqs.hostfree.html，就能看到交互式花冠图了。里面有专门的snapshot按钮用于截图。
在这里羡慕一下Jennifer Lu同学，能够在约翰斯霍普金斯大学进行宏基因组研究，并为我们带来了功能如此强大的软件和配套的工具。我也就只能在这所大学医学院的门口徘徊许久，喟叹而去，然后误入旁边的黑帮领地，惊恐万状，所幸全身而退。
除了使用seqs.hostfree.tsv来生成花冠图，seqs.hostfree.kra这个报告文件也可以，就是需要使用命令节选一下文件中的特定列：

# 同样需要注意命令的大小写。
$ ktImportTaxomomy -q 2 -t 3 seqs.hostfree.kra -o seqs.hostfree.html
# 或者是这样：
$ ktImportTaxonomy -m 3 -t 5 seqs.hostfree.kra -o seqs.hostfree.html

两种不同的方法不过是选用了报告数据中不同的列作为花冠图的参量，这个在GitHub论坛中都有用户进行过经验分享。如果我没记错的话，seqs.hostfree.kra的第2列和第5列应该是物种的编码，第3列是序列数。程序会自动对这些物种编码进行检索匹配，最终生成图像。然而为什么我很少用这种方法？因为这个文件里的物种都是用编码形式表示，这就意味着图像生成时需要使用数据库将其转换成物种名称。该命令运行时确实会下载并调用自己的数据库，但这也就意味着有些没有列入数据库编码的物种可能无法被识别出来。前一种方法直接使用文件当中的文字内容，虽然没有对这两种方法进行对比，但想想多少会全面一点儿吧。

序列组装

宏基因组目前普遍商用的方法还是采用二代测序，读长很短，需要将这些零碎的片段进行组装，才能一窥全貌。这里常用的还是两个大家熟悉的软件，MEGAHIT和SPAdes。

MEGAHIT序列组装

这个软件速度挺快，不过组装效果貌似稍微低一点，因此比较适合我这种配置比较低的NUC。使用方法也很简单：

$ megahit -r seqs.hostfree.fastq -o assemble_megahit

注意这里我使用的是-r选项导入数据，因为我是使用的单端序列。如果是双端就需要使用-1和-2选项了。组装之后的序列是在assemble_megahit目录下的final.contig.fa文件。

SPAdes序列组装

这个软件速度虽然相对比较慢，但是组装质量比较好。我的NUC也能比较快的跑完流程，估计是因为本身的数据量不够大，毕竟只用了处理之后hostfree的数据。使用方法也很简单：

# 这里使用-t选项调动了全部12个线程。
$ spades.py -s seqs.hostfree.fastq -o assemble_spades -t 12

同样注意这里使用的是-s来导入单端序列，双端也是使用-1和-2选项来导入。组装之后的序列是在assemble_spades目录下的contig.fasta文件。当然也同时生成框架序列的scaffolds.fasta文件，没怎么用过。

功能注释

功能注释是宏基因组分析里面要回答的第二个重要问题：在做什么。

HUMAnN功能注释

HUMAnN是一个非常强大的功能注释工具，我用它主要是觉得确实省事。一次性就搞定基因家族和相关通路，实在是不能太方便。也有很多其他的注释工具，这只是用得最顺手的而已。我在前一篇博文中提到了HUMAnN的安装，过程确实非常波折，希望我的分享能够帮助到更多的朋友。HUMAnN当中集成了MetaPhlAn这个物种注释工具，但是我使用过程当中觉得这个确实不如Kraken方便，而且不够稳定，数量也不够多，所以在此只是使用HUMAnN当中功能注释的部分。
这个软件的使用也是非常方便：

# 我已经使用humann_config将线程自动调整为全部12个，输出目录为工作目录下的humann3_out。
# HUMAnN3非常费时，主要是UniRef的数据库比上一代大了很多。可以使用--memory-use选项将所有的内存全部使用上应该能适当提速（默认是最小）。
$ humann -i seqs.fastq.gz -o humann3_out --memory-use maximum

注意这里是使用的CleanData数据seqs.fastq.gz进行分析。后续还有一些功能的挖掘，就需要使用组装之后的序列了。
如果是多个样本，就需要使用HUMAnH自带的一些工具进行整理，这个在以后关于后续处理的博文中再进行分享。

RGI耐药注释

细菌耐药是我在临床当中非常关心的问题，也是针对微生物相关的宏基因组分析当中的重要内容。有的公司会在报告中写这部分内容，但有的公司不会，或者不会常规分析，所以拿到数据我会自己比对一下。
耐药注释使用的是RGI工具，比对的是McMaster大学的CARD数据库。使用也很方便：

# 这里使用-n选项调用12个线程，导入的是之前megahit组装的序列。
$ rgi main -i assemble_megahit/final.contig.fa -o assemble_megahit/rgi -n 12

最后生成的是在assemble_megahit目录下的rgi.txt文件。没错，软件自动加文件名后缀，虽然是个txt文件，但实际上是一个tsv用Tab分隔的表格，可以用WPS或者Excel打开观看。同样的方法也可以处理assemble_spades目录下的contig.fasta文件。

VFDB毒力注释

细菌毒力也是我在临床工作中关心的问题。这里使用的是中国医学科学院病原生物学研究所维护的VFDB数据库。从网站上下载数据库，我选择完全版蛋白数据库。下载之后放在~/vfdb/目录下解压缩，里面有VFDB_setB_pro.fas和VFDB_setB_nt.fas两个序列文件，使用DIAMOND建立索引：

$ diamond makedb --in ~/vfdb/VFDB_setB_pro.fas --db ~/vfdb/vfdb_pro
$ diamond makedb --in ~/vfdb/VFDB_setB_nt.fas --db ~/vfdb/vfdb_nt

结束之后，就建立了vfdb_pro.dmnd和vfdb_nt.dmnd两个索引文件。现在进行比对。

$ diamond blastx --threads 12 --db ~/vfdb/vfdb_pro --query assemble_spades/contig.fasta --out vf.tsv

选项–out已经指定了输出文件是vf.tsv，结束之后就能打开看里面的比对结果。

后记

写完初稿已是深夜，窗外暴雨连绵。身在这样的一个开源互联的网络时代，我才能够有机会再次拿起键盘重启代码之旅。衷心感谢在网络上无私分享的每个人，正是这样的氛围，让我能够写下上面这些内容，加入到这样一个分享的行列当中。感谢这个时代。希望能够借此机会与更多的朋友进行交流，共同提高进步。文章当中存在的问题还希望大家及时指出，我会尽快修改完善。至于后续的数据处理，我会在以后单独写一个博文进行整理分享。
我爱着一个始终还存在的地球，一个乞丐的地球。它曾与极度的幸福结缘，又曾与极端的危险相依；既五光十色，又杂乱无章，它奇妙地存在着许多可能性。我一次又一次地战胜这个地球，由于它的美丽而疯狂，迷恋着它的模样，被强权所威胁，却又不被战胜。继续走吧，姑娘，向前走吧，孩子，我们交给了死亡，但是仍然活着。我的第二次获得的恩赐啊，你现在跟我一起迁徙。巴比伦，模糊而苍白，与它那石与钢造的、不断增高、抗拒着倒塌的塔楼一起崩溃；我们正急匆匆穿越风暴，后面骑兵在追踪，箭矢在射击。我们踩着沙子走，贴在斜坡上，脸被晒黑了。在我们前方，一个新国家远远地在朦胧中浮现，在银光中蒸腾。充满新的迫害，充满新的希望，充满新的歌唱！——《天使来到巴比伦》

宏基因组单个样本数据处理流程笔记相关推荐

青年生命科学论坛报告：扩增子和宏基因组数据分析与可视化流程—刘永鑫(北京210606)...
感谢中科院动物所青促会组织的第三届青年生命科学论坛的邀请,参加本次大会,并和微生物所王军老师共同负责了<微生物组>专题的召集工作.感谢11位微生物组专题报告人的辛苦准备和分享. 现将本次1 ...
宏基因组扩增子最新分析流程QIIME2：官方中文帮助文档
本网对Markdown排版支持较差,对格式不满意的用户请跳转至或"宏基因组"公众号阅读: 注:文为蓝色字均为文章链接,可点击直达写在前面 **声明:本文为QIIME2官方帮助文 ...
宏基因组扩增子2分析流程：中文首发,史上最详系,零基础自学
本网内容首发"宏基因组"公众号,更佳阅读体验.更多相关文章,欢迎点我跳转至公众号阅读写在前面之前发布的<扩增子图表解读>系列,相信关注过我的朋友大部分都看过了(链接 ...
宏基因组大数据分析的质量控制流程规范
宏基因组大数据分析的质量控制流程规范郑广勇1,杨桢1,曹瑞芳1,刘婉2,李亦学1,2,张国庆1,2 1. 中国科学院上海生命科学研究院生物医学大数据中心,上海 200031 2. 上海生物信息技术研 ...
你想要的宏基因组-微生物组知识全在这(2023.01)
欢迎点击上方蓝色"宏基因组"关注我们! 宏基因组/微生物组是当今世界科研最热门的研究领域之一,为加强宏基因组学技术和成果交流传播,推动全球华人微生物组领域发展,中科院青年科研人员创 ...
你想要的宏基因组-微生物组知识全在这(2022.12)
欢迎点击上方蓝色"宏基因组"关注我们! 宏基因组/微生物组是当今世界科研最热门的研究领域之一,为加强宏基因组学技术和成果交流传播,推动全球华人微生物组领域发展,中科院青年科研人员创 ...
你想要的宏基因组-微生物组知识全在这(2023.3)
欢迎点击上方蓝色"宏基因组"关注我们! 宏基因组/微生物组是当今世界科研最热门的研究领域之一,为加强宏基因组学技术和成果交流传播,推动全球华人微生物组领域发展,中科院青年科研人员创 ...
你想要的宏基因组-微生物组知识全在这(2023.7)
欢迎点击上方蓝色"宏基因组"关注我们! 宏基因组/微生物组是当今世界科研最热门的研究领域之一,为加强宏基因组学技术和成果交流传播,推动全球华人微生物组领域发展,中科院青年科研人员创 ...
你想要的宏基因组-微生物组知识全在这(2023.4)
欢迎点击上方蓝色"宏基因组"关注我们! 宏基因组/微生物组是当今世界科研最热门的研究领域之一,为加强宏基因组学技术和成果交流传播,推动全球华人微生物组领域发展,中科院青年科研人员创 ...

宏基因组单个样本数据处理流程笔记

宏基因组单个样本数据处理流程笔记

前言

数据预处理

质量控制

去除接头序列

去除宿主序列

物种注释

Kraken2注释

Krona制图

序列组装

MEGAHIT序列组装

SPAdes序列组装

功能注释

HUMAnN功能注释

RGI耐药注释

VFDB毒力注释

后记

宏基因组单个样本数据处理流程笔记相关推荐

最新文章

热门文章