高通量测序的数据处理与分析(二)-宏基因组4
基于read的物种注释
基于read的宏基因组物种注释的软件有很多，如 karen， kaiju , metaplhan。本文主要介绍如何使用 metaplhan进行基于read的物种注释。这种基于read的物种注释方法不依赖于contigs的组装质量，直接通过read和参考基因组比对得到丰度，从原理上讲会比基于bin模式的物种注释快上许多，但是这种基于read的比对模式高度依赖参考基因组的准确性，完整性。适用于人类相关微生物，小鼠相关微生物或者其他研究较为深入的微生物环境。而比较复杂的环境微生物可能用这种方式注释出来的物种会不太准确。
安装metaphlan
# 创建一个隔离的conda环境
conda create --name metaphlan python=3.8
# 安装metaphlan
conda install -c bioconda metaphlan
pip install metaphlan # 如果conda安装不上就用pip
mamba install -c conda-forge -c bioconda metaphlan # 如果还安不上用mamba
# 安装数据库, 注意数据库无需重复安装，在使用时用 --bowtie2db参数进行指定就行
metaphlan --install --bowtie2db [path/of/databasefile]
若无法使用命令自动下载，则手动下载 .tar、.md5 和 mpa_latest 文件放到自定义的数据框位置即可。数据库下载地址：Segatalab FTP
metaphlan的使用
metaphlan [ERR_num_1.fastq],[ERR_num_2.fastq] --nproc 5 --input_type fastq --bowtie2db [path/of/databasefile] --bowtie2out metagenome.bowtie2.bz2 -o [path/of/out/ERR_num.txt] -t rel_ab_w_read_stats
--nporc ： 指定核心数
--input_type ：指定输入文件类型
--bowtie2db ：指定数据库位置
--bowtie2out ：保存中间结果，可以提速
-t:
rel_ab # 输出物种的相对丰度，
rel_ab_w_read_stats # 同时输出物种的相对丰度和每千碱基的read数，这个read数可以直接用于DESeq2求差异
reads_map # 仅输出read
marker_counts # 非标准化marker counts，输出的SGB编号的marker counts，什么是marker counts不太清楚
....
默认是rel_ab
metaphlan的结果
#mpa_vOct22_CHOCOPhlAnSGB_202212
#/home/baizy2018/anaconda3/envs/metaphlan/bin/metaphlan ERR2759385_1.fastq,ERR2759385_2.fastq --nproc 12 --input_type fastq --bowtie2db /home/baizy2018/wxl/database/metaphlan --bowtie2out ./metagenome.bowtie2.bz2 -o ./ERR2759385.txt -t rel_ab_w_read_stats
#2653050 reads processed
#SampleID Metaphlan_Analysis
#estimated_reads_mapped_to_known_clades:2213524
#clade_name clade_taxid relative_abundance coverage estimated_number_of_reads_from_the_clade
k__Bacteria 2 100.0 0.89774 2213524
k__Bacteria|p__Actinobacteria 2|201174 75.36809 0.67661 1573892
k__Bacteria|p__Firmicutes 2|1239 18.14211 0.16287 432143
k__Bacteria|p__Bacteroidetes 2|976 5.27663 0.04737 180609
k__Bacteria|p__Proteobacteria 2|1224 1.19424 0.01072 26406
k__Bacteria|p__Verrucomicrobia 2|74201 0.01894 0.00017 474
....
其中
- clade_name： 分类名，这些clade_name有分到
t__级别，t级别是株水平的分类[参考 1]。 - clade_taxid： 分类ID
- relative_abundance： 物种相对丰度
- coverage： 覆盖度
- estimated_number_of_reads_from_the_clade： counts数
有些结果可能会有 additional_species 列，这些列是表示次选的物种[参考 2]