宏基因组宏基因组

2018-08-20 Metagenome分析教程

2018-08-20  本文已影响0人  电泳菌要努力
pic

前言

之前看panda姐的blog时，发现介绍了篇关于Metagenome的分析教程，刚好就想着试着翻译下，作为我的第一篇文章，顺带练练我MarkDown语法的训炼吧。

Introduction

目前微生物学在宏基因组方面最大的挑战是需要将宏基因组学和序列分析纳入传统微生物学的训练中。Sean Eddy(Howard Hughes医学中心的计算机生物学家）非常6的对生物学上的高通量测序和其训练在他的keynote进行了讲解。

为了方便微生物学家组装宏基因组，我们这次补充了如何评估宏基因组里参考序列（如，gene,contigs,等等）丰度的指导。我们的方法包括以下所有的参考1.可获得的基因组参考 2.由宏基因组组装来的。通常来说，你想完成该指导，或者大多数的宏基因组组装，需准备以下东西：

1.可用的服务器。大多数宏基因组组装需要的内存是大多数研究人员的个人电脑达不到的。在本次讲解中，我们建议您使用公开的Amazon EC2样本服务器，所有人经过登记后都可以注册。

2.一个宏基因组数据库。我们本次讲解选择使用WGS数据库的HMP MOCK种群，鉴于其的可使用性，合适的大小，以及参考基因的可获得性。该数据库是22种种类，其基因被从自然界分离，合并并测序，而来的模拟宏基因组。

3.负责组装，read mapping,和基因注释的软件，我们将会展示将这些软件安装在Ubuntu系统的服务器上。

0.和我们登陆同一个页面

本次教程的第一步是提供了一面向所有用户，无论你使用寿命电脑的可使用以下全部命令的服务器。为了达到这个目的，我们需要使用云计算。具体来说，亚马逊网络服务弹性计算云。租用该网络服务，你需要申请一个信用卡，并使用该信用卡付费。该网络服务收费很合理。你完成该教程大约在四小时以内，花费低于1美元。（译者注：国内可以使用阿里云或百度云完成该教程，注意，在校学生如果没有接触过Linux操作系统，可以先行购买10元/月的学生优惠对Linux命令与软件安装进行学习，但是该服务器跑RNA-seq分析都跑不动啊。）

当你注册好了亚马逊网络服务器，你需要按照网站说明启动云“案例”或服务器。作为本次教学，我们建议你使用Data Science Toolbox。在运行前有几点需要注意（不用Amazon的话，以下意义不大）：

1.选择"in the cloud"的运行指导。

2.你可以使用任一AMI（系统映像），但我们建议你使用US EAST,ami-d1737bb8。（这个国内的看看就好了。)

  1. 不要忘记Data Science Toolbox指导里的第二条规制：添加“自定义TCP规制”里端口为“8888”，源为“Anywhere"。

4.从头至尾的完成到了step4后。当你到达step5时，参考下文。

如果你用的是Mac或Linux操作系统登录远程服务器出现问题：

。检查你的密钥权限是否被更改（就是那个以*.pem结尾的那个）

。确定你在跑ssh命令远程登录服务器和你的安全文件在同一个文件夹下。

当你可以成功的跑类似于以下命令时，你便可以成功登录服务器（其中的名字和EC2地址是用你自己的）

    $ ssh -i MyKeyPair.pem ubuntu@ec2-XX-XX-XX-XXX.compute-1.amazonaws.com

此时你的命令行变为类似于如下所示：

ubantu@ip-10-181-106-120:

在运行本教程前你还需要做几件小事，

逐条复制粘贴以下命令在你的命令行中，并每条后用ENTER运行：

    cd /mnt

    sudo git clone https://github.com/germs-lab/frontiers-review-2015.git

之后输入下列命令行，并在提示出现后输入笔记本密码：

    dst setup base

之后再复制粘贴以下命令：

    sudo ipython notebook --profile=dst --notebook-dir=/mnt/frontiers-review-2015

这样会为本次教程启动一个Ipython Notebook。离开终端，打开你的互联网浏览器，最好是Google Chrome。你还需在上面登录你的EC2实例的公共DNS地址，如"ec2-XX-XX-XX-XXX"。如果你不知道你的DNS地址，你可以随时在你的AWS EC2控制板上查看。

在你的浏览器上，导航至https://ec2-XX-XX-XX-XXX:8888。几乎所有浏览器都会提醒你这是个不安全的地址，不需要紧张。在Chrome浏览器上，你可以点击“高级”选项，然后再点击“继续进程”。然后输入密码（之前你输入的笔记本密码），然后你就愉快的看见一个包含着一个叫”frontiers-nb-2015“文件的笔记本了。

1.如何使用这个IPython笔记本

IPython笔记本在共同训练生物信息学时很有用。这些笔记最近被Nature News ( http://www.nature.com/news/interactive-notebooks-sharing-the-code-1.16261 和 http://www.nature.com/news/programming-pick-up-python-1.16833)报道。

2.下载本次教程的数据

我们从NCBI Short Read Archives(SRA)下载HMP模拟宏基因组作为本次教程的开端。NCBI的SRA中存储了许多公开的宏基因组。使用sratoolkit软件是最简单的获得SRA数据的方法。如果你知道你需要的SRA run ID（比如SRR172903)，你可以用专门的程序下载该数据并转为”fasta“或者”fastq"格式序列。

    wget http://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.4.5-2/sratoolkit.2.4.5-2-ubuntu64.tar.gz

    tar -xvf sratoolkit.2.4.5-2ubuntu64.tar.gz

你现在使用ls命令可以发现多出来个包含软件的文件夹。你同时在我们工作的notebook上看见该文件夹。

    ls

现在，我们可以使用安装好的sratoolkit程序去下载HMP模拟数据库的"fastq"格式了。（大概会花费1到2分钟。耐心在宏基因组的研究中很有必要，好消息就是你依靠“租用”的云服务器工作，不会占用你个人电脑内存的运算能力--这样等待的过程中，你可以做自己的事。你会注意到有个“Kernel busy"会显示在"登出"按钮下方屏幕的右上角处。

    sratoolkit.2.4.5-2-ubuntu64/bin/fastq-dump SRR172903

3.质量控制

有很多种方法可以去测定你组装的序列数据的质量。首先，你可以查看你测序reads的质量分数，如果必要，你可以使用分数不满足的reads，之后进行修剪。有大量的软件可以进行修剪，一些软件还包含着很好的教程，包括FastX Toolkit(http://hannonlab.cshl.edu/fastx_toolkit/ 和 http://khmer-protocols.readthedocs.org/en/v0.8-1/metagenomics/1-quality.html)，FastQC(http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ 和 http://ged.msu.edu/angus/tutorials-2013/short-read-quality-evaluation.html)与及Sickle(http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ 和 http://ged.msu.edu/angus/tutorials-2013/short-read-quality-evaluation.html)。

你之前下载的“fastq”格式的序列文件，其每个序列read由四行展示，我们做个快速的浏览：

    head -n 4 SRR172903.fastq

. 第一行（开头是“@SRR172903.1”)是read的标识符，它通常是展示read的ID，测序仪器的一些关于如何获取序列的信息。

. 第二行是DNA序列。

. 第三行的内容以第一行相同，只是把"@"换成了"+"，有时在部分数据库内只有一个"+"。

. 第四行是DNA序列每个碱基对的质量评分信息，注意该部分是与DNA序列长度一致的，且该质量得分是基于ASCII码字符分数（其具有测序技术确定偏移，Illumina目前的偏移是64，例如，ASCII编码 64 = 0 Phred 得分）。质量得分等于 -10*log(p), p是该碱基错误的概率（比如，如果Q=20, p=0.01,就是有1%的几率这个碱基是错误的。

这次教程中，我们会将超过序列内超过50%的碱基Phred得分低于33分的read移除。移除过程将用Fastx-Toolkit完成，该软件可以进行多种质量控制（比如adapter剪切）。我们首先来下载，解压和安装这个软件吧。

    wget https://github.com/agordon/fastx_toolkit/releases/download/0.0.14/fastx_toolkit-0.0.14.tar.bz2

    wget https://github.com/agordon/libgtextutils/releases/download/0.7/libgtextutils-0.7.tar.gz

    tar -xvf fastx_toolkit-0.0.14.tar.bz2

    tar -xvf libgtextutils-0.7.tar.gz

    bash fastx_install.sh

现在，我们可以开始质量过滤了，并将过滤后的文件储存为SRR172903.qc.fastq。

FASTQ质量控制

$ fastq_quality_filter -h usage: fastq_quality_filter [-h] [-v] [-q N] [-p N] [-z] [-i INFILE] [-o OUTFILE]

version 0.0.6

[-h] = This helpful help screen.

[-q N] = Minimum quality score to keep.

[-p N] = Minimum percent of bases that must have [-q] quality.

[-z] = Compress output with GZIP.

[-i INFILE] = FASTA/Q input file. default is STDIN.

[-o OUTFILE] = FASTA/Q output file. default is STDOUT.

[-v] = Verbose - report number of sequences.

If [-o] is specified, report will be printed to STDOUT.

If [-o] is not specified (and output goes to STDOUT),

report will be printed to STDERR.

    fastq_quality_filter -q 33 -p 50 -i SRR172903.fastq > SRR172903.qc.fastq

4.多样化检查---“都有谁在”的分布情况

宏基因组测序的优势在于其无需微生物培养就能够量化微生物在环境中分布的能力。通常大多数研究需要了解其分类多样性（尤其是对16s rRNA基因测序的研究中）。多样性也能在宏基因组的特定序列模式的表示下被测定。比如，可以量化宏基因组中独特的核苷酸k-mers。这些k-mers也能用作组装宏基因组，重叠的k-mers被用于指出read应该被连接在一起。k-mers的多样性可以给你用来研究你样本的多样性。同时，由于组装对比中每一个k-mer需要与所有的k-mers进行对比，所以大量的k-mers存在需要更多的计算机内存。Miller等人对k-mers和组装做了个很好的综述。

（需要注意的是，16S rRNA扩增子测序是一种靶向的方法，在本教程中默认其不属于宏基因组测序。鸟枪法宏基因组测序是将种群中全部细胞的DNA提取出来进行测序。而靶向测序扩增特点的基因座并独立测序。Sharpton等人在宏基因组分析上有很好的综述。）

现在我们需要做的是安装khmer--它包括了一套khmer和pre-assembly工具。我们在这里用于k-mer的计数。一旦你开始运行一下脚本，你能够开始使用khmer的许多工具了。

    ls

    bash khmer-install.sh

接下来的脚本包含在khmer包中，并可以估算数据集的k-mers的唯一总数。用该数据可以a.)明确宏基因组组成的多样性，例如用于组装的细菌基因组，b.)比较复数宏基因组之间的k-mer多样性，c.)探讨k长度选择对组装的影响。

然后为了评估不同k值（17，21，25，29，33，37）的k-mers的唯一总数，运行下面的脚本。该脚本输出唯一k-mers但也将将其保存在unique_count文件中。（该过程在一个大型服务器实例中大概耗时15分钟，在超大型实例中耗时8-10分钟。）

    python unique-kmers.py -R unique_count -k 17 SRR172903.qc.fastq

    python unique-kmers.py -R unique_count -k 21 SRR172903.qc.fastq

    python unique-kmers.py -R unique_count -k 25 SRR172903.qc.fastq

    python unique-kmers.py -R unique_count -k 29 SRR172903.qc.fastq

    python unique-kmers.py -R unique_count -k 33 SRR172903.qc.fastq

    python unique-kmers.py -R unique_count -k 37 SRR172903.qc.fastq

现在你可以看见这个文件，第一栏显示k-mer长度，第二栏显示宏基因组中所要求的k长度的数量。如果你有复数的基因组，你可以通过比如数据库里的k-mers总数比较其多样性。你可以使用命令“cat”来查看文件里的结果。

    cat unique_count

5.得到基因的覆盖描述：我的宏基因组目前有哪些基因？

大多数的宏基因组分析需要评估参考基因的丰都（比如，来自基因组或者自己组装的宏基因组）。本教程将覆盖参考基因可以得到或不可以得到（需要de nove组装）的两种情况。

6.例子1--有可获得的参考基因

对于模拟的HMP宏基因组，HMP已经对用于模拟数据库的被分离的基因组进行过测序。这些基因组的列表可以在HMP网站上获得，同时我们将他上传到了github仓库---一个用来协作分享数据和代码的工具。本次教程中用以下命令就可以下载该数据。

    cat ncbi_acc.txt

下面的命令用于下载上面那个列表中的每一个ID的基因组进入一个叫“genomes"的字典。

    python fetch-genomes-fasta.py ncbi_acc.txt genomes

7.估计组装contigs的丰度

为了评估你基因组里参考基因和基因组的表现，你可以用mapping软件将reads比对到参考基因上（比如，Bowtie2,BWA等等）。在本次教程中，我们将在我们的服务器上安装并使用Bowtie2。我们之后会mapping我们的宏基因组到单参考基因组上（我们之前下下来的。）。

    wget http://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.5/bowtie2-2.2.5-linux-x86_64.zip

    unzip bowtie2-2.2.5-linux-x86_64.zip

我写一个的脚本能够自动将一系列read map到给与的参考基因上并输出一个包含能够被map上的reads数量的文件。为了使用该脚本，我们需要安装samtools。samtools帮我们与一种超级压缩文件samfile交互，该文件可以有效的从mapper中储存mapped信息。

    apt-get install samtools

我们提供一个简单的操作过程来将reads map到参考基因组上。该过程操作如下述步骤：

.制作你的参考基因索引

.将你的reads map到你的参考基因索引中（用默认的bowtie参数）

.用Samtools评估你的read的map数，没有map上的数，并且提供一个制表符分割文件，每行一个参考基因名和序列长度，map上的read和未map上的read。

该过程大概花8到10分钟。

    bash bowtie.sh genomes/NC_000913.2.fa SRR172903.qc.fastq

现在我们可以看看基因组NC_000913上map和没有map到的read了。我们同时有了个文件可以展示参考基因名（第一列），参考基因长度（第二列），map上的read(第三列）和未map上的read(最后一列）。你可以点击这里查询关于samtools的更多信息。

    cat reads-mapped.count.txt

    cat reads-unmapped.count.txt

如果你想要挑战一下更难的任务，你可以试试map这个宏基因组到基因组文件夹提供的所有参考基因组上。试试串联所有基因组到一个文件夹里吧（命令是"cat genomes/*fa >> all-genomes.fa"）并运行对所有genomes.fa运行脚本，而不仅仅是对NC_000913.2.fa。

8.例子2--De novo组装参考基因

HMP模拟宏基因组的组装

组装是将宏基因组里有望来自同一基因的reads的重复部分合并，使其成为一个较长的，连续的序列（通常称之为一个contig）。其优点是提供一个更长的序列，使其能过在未来作为参考基因（之前未知的），减少分析的数据量，并提供一个不依赖过去知识的参考。

选择使用什么组装方式并不容易且存在争议。请记住重要的一点，组装是对你的数据的公认的假设表达。组装本身只是第一步，接下来需要对其精确性和实用性进行评估。对于大多数的组装情况，需要将测序的reads和参数输入组装软件。本次教程我们将使用2014年发布的Megahit(Li et al.,2015,https://github.com/voutcn/megahit)进行组装。Sharpton的综述(Sharpton, 2014)同样也介绍了许多很好的宏基因组组装流程和方法。

为了减少所需的内存，常常需要事前统一宏基因组中k-mers的分布。删除无需组装的额外信息和可能错误的reads可以优化组装(http://arxiv.org/abs/1203.4802)。这些脚本与教程你可以从 http://ged.msu.edu/angus/diginorm-2012/tutorial.html 中获得。

本次教程中，我们使用Megahit组装宏基因组，所以，让我们先安装megahit。

    bash install-megahit.sh

该过程大约15分钟，然后我们会获得一个叫做”megahit_assembly"的文件夹。你可以用在这里阅读其参数，比如 --memory是指定其可以使用的最大内存。

    megahit/megahit --memory 10e9 -l 250 --k-max 81 -r SRR172903.qc.fastq --cpu-only -o megahit_assembly

为了看看这个组装，让我们在最终contigs所在的megahit_assembly.final.contigs.fa文件下运行khmer组装总程序。我们统计下contigs大于等于200bp的。

    python khmer/sandbox/assemstats3.py 200 megahit_assembly/final.contigs.fa

9.评估contigs丰都

当组装结束是，你有一组可以评估宏基因组丰都的参考contigs。评估方法与之前使用参考基因组的方法相同。

这些会花20分钟。

    bash bowtie.sh megahit_assembly/final.contigs.fa SRR172903.qc.fastq

你可以如同将reads map到NCBI基因组一样的查看下列map结果。

    cat reads-mapped.count.txt

    cat reads-unmapped.count.txt

    cat reads.by.contigs.txt

10.注释组装后的contigs

测序通常是来确定你样品里的“who"和"what"。在我们的例子里，我们知道HMP模拟种群来源的一组基因组（我们在之前下载的。）。将未知序列与已知参考基因进行对比的是最为流行的工具中的其中之一--The Basic Local Alignment Search Tool（BLAST)。为了明确我们的ontigs的来源，我们将组装的contigs和HMP模拟种群的基因组进行比对。

我们首先下载BLAST软件。鉴于序列数据库容量的不断增加，也可以考虑使用目前在注释中更有效率的新工具，比如Diamond(https://github.com/bbuchfink/diamond/, http://dx.doi.org/10.1038/nmeth.3176)。

    wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.2.30/ncbi-blast-2.2.30+-x64-linux.tar.gz

    tar -xvf ncbi-blast-2.2.30+-x64-linux.tar.gz

现在我们可以为BLAST做一个可搜索的数据库了。首先，我们要把在基因组字典里的所有基因组连接到同一个文件里。

    cat genomes/*fa >> all-genomes.fa

    ncbi-blast-2.2.30+/bin/makeblastdb -in all-genomes.fa -dbtype nucl -out all-genomes

    ncbi-blast-2.2.30+/bin/blastn -db all-genomes -query megahit_assembly/final.contigs.fa -outfmt 6 -out contigs.x.all-genomes.blastnout

上面的命令比对了每一个query（每一个在组装好的final.contings.fa文件里的序列）与每一个序列（比如all-genomes.fa中的基因组）。-outfmt告诉该程序将结果以制表符格式存在 -out 的 contigs.x.all-genomes.blastnout 文件中。

让我们瞧一眼这个文件的前10行。你会看见query(contig)和hit(genome)后的同一性百分比，对比长度，错配计数，间隙开放计数，query开始位置，query结束位置，subject开始位置，subject结束位置，E-值和比特得分。

    head -n 10 contigs.x.all-genomes.blastnout

根据你的研究问题，对开放阅读框（ORFs)的注释可能比对contigs序列的注释更有意义。本次例子中，存在多个可以使用的ORF caller((e.g., FragGeneScan, http://nar.oxfordjournals.org/content/early/2010/08/29/nar.gkq747.abstract 和 Metagene, http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1636498/)。我们可以用FragGeneScan从我们的contigs调用ORFs。首先依旧是按如下方式下载，安装软件，然后从我们的contigs中调用ORFs:

    wget http://downloads.sourceforge.net/project/fraggenescan/FragGeneScan1.19.tar.gz

    tar -xvf FragGeneScan1.19.tar.gz

    bash fraggenescan-install.sh

我们将在组装的contigs上运行FragGeneScan，假设它符合“完整”基因组序列的训练特征（在他们的文档中，这等于完整的基因组序列或短序列读数而没有测序错误）。

    FragGeneScan1.19/FragGeneScan -s megahit_assembly/final.contigs.fa -o final.contigs.orfs.fa -w 1 -t complete

ORFs调用后最终组成为一个FASTA文件，存在final.contigs.orfs.fa.faa（氨基酸）和final.contigs.orfs.fa.ffn（碱基）中，你可以根据以上contigs描述的数据对你选择的数据库进行描述。

11.以后的事

现在你知道获得以下信息的所有该知道的事了：

1，序列丰都信息：序列（比如contig)和丰都（比如，可map的read数）

2，序列注释信息：序列（比如contig)和NCBI基因组

你会注意到这与16S rRNA扩增子分析相似，你可以在其中获得OTU丰度表和OTU最佳命中注释。对于宏基因组分析，得到这些信息你可以用于进一步的分析和可视化软件包，如R语言中的PhyloSeq包。

译者最后说下

本文是我在学习宏基因组的过程中进行翻译的，所以有很多地方还有错误，同时你们也会发现包括read,contig,map之类的词...我完全不知道它们的中文是什么呀，请大家原谅，如果发现有什么地方不对多多指教，在原文中，每一个linux脚本代码前都有”！“，我学了那么久的linux，没见过这样的用法，我一开始担心是格式转换问题，所以我默默的删去了，最后原文有写为什么加”！“，好像是Ipython notebook的内置语法，我没太搞清楚。这也是我blog的第一篇文章，以后看到关于肠道微生物，宏基因组，Python和R语言的好文章，或者有写什么总结，会持续发在我的Git上的，emmmm....虽然我想看得人不会很多，但是还是写一下吧，如果有人看当然更有动力啦。啊哈哈哈哈哈~

上一篇下一篇

猜你喜欢

热点阅读