免组装宏基因组群落分析
宏基因组数据可以不经组装，直接将测序获得的reads比对到公共数据库中，利用比对到的数据库序列的物种归属信息对reads进行物种分类，从而快速获得群落的物种组成信息。
宏基因组分析Pipeline
免组装宏基因组群落分析
更新中……
01
KAIJU
KAIJU（http://kaiju.binf.ku.dk/）是一个对宏基因组高通量测序数据进行物种分类的工具，Kaiju将所有的reads翻译成氨基酸序列，然后在蛋白质数据库（NCBI RefSeq或者NR）中搜寻这些序列，来发现最精准的比对（maximum exact matches，MEMs）。
在使用前，需要使用参考蛋白质数据库构建Kaiju的数据库Index，首先从NCBI的FTP下载参考基因组和taxonomy文件，可以使用Kaiju的makeDB.sh命令自动下载并构建Index，如下所示：
mkdir kaijudb
cd kaijudb
nohup ../bin/makeDB.sh -e &
#具体参数如下所示：
-r：只下载完整的细菌和古菌基因组数据从NCBI RefSeq
-p：下载proGenomes数据库（http://progenomes.embl.de/）所有蛋白序列
-n：下载NCBI NR数据库并提取属于细菌、古菌和病毒的蛋白质序列
-e：下载NCBI NR数据库并提取属于细菌、古菌、病毒以及真菌、真核微生物的蛋白质序列
-m：下载海洋宏基因组网站（Marine Metagenomics Portal，MMP，https://mmp.sfb.uit.no/）海洋参考数据库MarRef和MarDB
-v：是-r与-p的补充，下载病毒的基因组数据
-t：构建Kaiju的Index时使用的核数，默认为5
--noDL：不下载数据库文件，使用已经存在的文件构建Index，例如-e参数会下载nr.gz与taxdump.tar.gz，如果存在已经下载的这两个本地文件，就可以不再重复下载Kaiju主程序使用方法如下所示：
$PATH/kaiju/bin/kaiju -t nodes.dmp -f kaiju_db.fmi -i reads.fastq [-j reads2.fastq]
参数含义如下所示：
-t：数据库nodes.dmp文件名，在kaijudb，需要加绝对路径
-f：数据库(.fmi)文件名，在kaijudb，需要加绝对路径
-i：reads文件名，fasta或者fastq格式
-j：双末端测序中第二个reads文件名
-o：输出结果文件名，如果没有设置，默认输出到标准输出
-z：程序运行的核数，默认为1
-a：使用最精准比对算法"mem"还是贪婪算法"greedy"，MEM算法只考虑精确匹配，不允许有错配，默认为greedy
-e：贪婪算法中最大允许错配数，默认为3
-m：最小比对长度，默认为11
-s：在贪婪算法中最低匹配得分，默认为65
-E：在贪婪算法中比对的最低E-value，与blastp类似
-p：输入序列为蛋白质序列
-v：输出备份信息下面对质控后的clean reads进行物种注释：
nohup kaiju -z 20 -t nodes.dmp -f kaiju_db_nr_euk.fmi -i screen.clean_1.fq -j screen.clean_2.fq -o clean.reads.kaiju.nr.out &结果中包含每条序列的taxid，但是并不包含物种名称，可以为注释结果添加各层级的物种名称：
nohup addTaxonNames -i clean.reads.kaiju.nr.out -t nodes.dmp -n names.dmp -o clean.reads.kaiju.nr.names.out -p -u &其中-u表示无法分类的reads不包含在输出结果中，-p表示输出完整的分类路径，也即不同层级的分类单元。
基于比对结果，可以在不同分类水平对reads进行binning，从而获得不同分类水平的群落结构信息，下面以门水平为例：
kaijuReport -i clean.reads.kaiju.nr.out -t nodes.dmp -n names.dmp -o clean.reads.kaiju.out.phylum.summary -r phylum
其中-r为分类水平，binning结果如下所示：
Krona（https://github.com/marbl/Krona/wiki）是一个很好的分层数据探索工具，通过可缩放的、多层的扇形图进行展示数据结构，krona数据可以通过KronaTools绘制图形，首先使用Kaiju的kaiju2krona工具生成krona数据：
kaiju2krona -i clean.reads.kaiju.nr.out -t nodes.dmp -n names.dmp -o clean.reads.kaiju.nr.krona然后利用KronaTools绘制图形：
ktImportText -o clean.reads.kaiju.nr.krona.html clean.reads.kaiju.nr.krona结果如下所示：
01
MetaPhlAn
MetaPhlAn（http://huttenhower.sph.harvard.edu/metaphlan/）是分析微生物群落物种（细菌、古菌、真核生物和病毒）组成的Pipline，它在宏基因组研究中非常有用，只需一条命令即可获得微生物的物种丰度信息，目前最新的版本是MetaPhlAn2（http://huttenhower.sph.harvard.edu/metaphlan2）。
MetaPhIAn的基本处理思想为首先将已知数据库的序列信息进行分析，最终形成每个物种独特的marker，然后将测序数据跟marker进行blast比对，确定物种类别，最终根据每个物种比对上的reads数目以及marker长度计算得到丰度。进化分支特异的maker（Clade-specific markers）需要满足在该分支内的基因组中是保守的，并且与分支外的基因组序列不相似。在MetaPhlAn中，物种分类准确性在于物种的基因组数据是否足够丰富，越丰富，marker的信息越准确，此处用2887个基因组数据进行的marker计算。计算速度的快慢取决于marker和reads的比对速度，此处用BLAST进行的操作。为此在MetaPhlAn2中，数据库进行了更新，使用了约17000参考基因组数据，(其中细菌和古菌13500、病毒3500、真核微生物，并用bowtie2进行reads和marker的比对。
软件下载解压后即可使用，主程序脚本为metaphlan2.py，第一次使用，程序会自动从https://bitbucket.org/biobakery/metaphlan2/downloads/下载数据库文件（mpa_v20_m200.tar）至安装目录的metaphlan_databases文件夹中，并进行校验、解压。软件使用方法如下所示：
metaphlan2.py [Options] INPUT_FILE OUTPUT_FILE
Options：
--input_type：输入文件格式，可选fastq、fasta、multifasta、multifastq、bowtie2out、sam，默认为automatic，也即自动检测
--bowtie2db：BowTie2数据库，默认安装目录下metaphlan_databases文件夹
--bowtie2_exe：BowTie2软件脚本的名称及路径，当其不在默认环境变量时可以使用此参数来指定
--bowtie2_build：bowtie2-build程序路径，当其不在默认环境变量时可以使用此参数来指定
--bowtie2out：储存BowTie2比对结果的文件名称，类似于bed文件
--tmp_dir：临时文件储存路径，默认为系统tmp dir
--tax_lev：相对丰度输出结果的分类层级，可选a（all taxonomic levels）、
k（kingdoms）、p（phyla only）、c（classes only）、o（orders only）、f（families only）、g（genera only）、s（species only），默认为a
--min_alignment_len：reads最小比对长度
--ignore_viruses：忽略病毒
--ignore_eukaryotes：忽略真核生物
--ignore_bacteria：忽略细菌
--ignore_archaea：忽略古菌
--sample_id_key：指定样品ID key，默认为#SampleID
--sample_id：指定sample ID，默认为Metaphlan2_Analysis
-s, --samout：输出的sam文件名称
--mdelim, --metadata_delimiter_char：不同分类层级之间的分隔符，默认为pipe，也即：k__Bacteria|p__Proteobacteria
--nproc：程序运行所使用的核数，默认为4
--stat_q：用于截断或缩尾统计的分位数
--stat：将markers丰度转换为系统发育分支丰度的统计方法，有以下几种（默认为tavg_g）：
avg_g：全部marker丰度的均值
avg_l：长度标准化的marker丰度的均值
tavg_g：总体marker丰度的截短均值，基于分位数--stat_q
tavg_l：长度标准化的marker丰度的截断均值，基于分位数--stat_q
wavg_g：总体marker丰度的缩尾均值，基于分位数--stat_q
wavg_l：长度标准化的marker丰度的缩尾均值，基于分位数--stat_q
med：长度标准化的marker丰度的中位数
INPUT_FILE：输入文件名称
OUTPUT_FILE：输出文件名称所谓截断均值就是去掉该分类层级最低的一部分taxonomy，使用剩余高丰度部分来计算相对丰度均值，而缩尾均值则是使用次低的部分来代替最低的部分。接下来使用metaphlan2对宏基因组clean reads进行分析：
nohup metaphlan2.py --nproc 20 --stat tavg_l --bowtie2out meta.bowtie2.bz2 --input_type fastq <(cat screen.clean_1.fq screen.clean_2.fq) > metagenome.txt &结果如下所示：
在metaphlan2安装目录下，utils文件夹中所包含的脚本merge_metaphlan_tables.py可以将不同样品的物种谱融合在一起，方便后续的比较分析，多个文件空格分隔，或使用通配符，如下所示：
merge_metaphlan_tables.py 1061_metagenome.txt 1171_metagenome.txt 1194_metagenome.txt 1237_metagenome.txt | sed 's/_metagenome//g' > merged_metagenome.txt结果如下所示：
此外脚本metaphlan_hclust_heatmap.py可以绘制热图来比较不同样本的物种构成差别，使用方法如下所示：
metaphlan_hclust_heatmap.py [options] --in INPUT_FILE --out OUTPUT_FILE
--in：使用utils/merge_metaphlan_tables.py获得的结果
--out：输出结果图像，后缀可选png、pdf、svg
-m：层次聚类方法，可选single、complete、average、weighted、centroid、median、ward，默认为average
-d：样品距离计算方法，用于样品聚类，可选euclidean、minkowski、cityblock、seuclidean、sqeuclidean、cosine、correlation、hamming、jaccard、chebyshev、canberra、braycurtis、mahalanobis、yule、matching、dice、kulsinski、rogerstanimoto、russellrao、sokalmichener、sokalsneath、wminkowski、ward，默认为braycurtis
-f：物种距离计算方法，用于物种聚类，可选euclidean、minkowski、cityblock、seuclidean、sqeuclidean、cosine、correlation、hamming、jaccard、chebyshev、canberra、braycurtis、mahalanobis、yule、matching、dice、kulsinski、rogerstanimoto、russellrao、sokalmichener、sokalsneath、wminkowski、ward，默认为correlation
-s：数据范围缩放方法，可选log
-x：热图一个小单元的宽度，对于非常大的图可设置此值
-y：热图一个小单元的高度，对于非常大的图可设置此值
--minv：热图展示的最小值（相对丰度），默认为0也即全部展示
--maxv：热图展示的最大值（相对丰度），默认为出现的最大值，设置成100则展示所有物种
--tax_lev：展示的分类层级，可选a（all taxonomic，混杂在一起，不推荐）、k（kingdoms）、p（phyla）、c（classes）、o（orders）、f（families）g（genera）、s（species），默认为s
--perc：百分数，只显示高丰度物种，默认为90
--top：数字，只展示排名前几的高丰度物种ordering based on --perc)
--sdend_h：样品树的高度，默认为0.1
--fdend_w：物种树的宽度，默认为0.1
--font_size：标签字体大小，默认为7
--clust_line_w：系统发育树线的宽度
-c：绘制热图的颜色设置：Accent、Blues、BrBG、BuGn、BuPu、Dark2、GnBu、Greens、Greys、OrRd、Oranges、PRGn、Paired、Pastel1、Pastel2、PiYG、PuBu、PuBuGn、PuOr、PuRd、Purples、RdBu、RdGy、RdPu、RdYlBu、RdYlGn、Reds、Set1、Set2、Set3、Spectral、YlGn、YlGnBu、YlOrBr、YlOrRd、afmhot、autumn、binary、bone、brg、bwr、cool、copper、gist_earth、gist_ncar、gist_rainbow、gist_stern、gray、hot、hsv、jet、pink、seismic、spectral、spring、summer、terrain、winter、bbcyr、bbcry，默认为jet下面对4个样品的融合结果进行作图，展示属水平所有物种的群落：
metaphlan_hclust_heatmap.py --tax_lev g -c Blues --in merged_metagenome.txt --out merged_metagenome.pdf作图结果如下所示：
—END —
腾讯云开发者
