一、metaphlan简介
MetaPhlAn 是用于二代测序物种分类的工具,只需一条命令即可得到宏基因组物种分类的列表,以及相对丰度信息,使用起来非常方便。目前最新版本是 metaphlan3.0。
metaphlan 是利用快速比对工具 bowtie2 与 Marker 基因集进行比对,而非全基因组序列,因此运行速度非常快。直接输入 kneaddata 质控过滤后的数据,即可快速获得微生物群体中种水平精度的组成,包括细菌、古菌、真核生物和病毒。
metaphlan 输出的文件格式是一种非常通用的列表格式,非常方便下游软件进行进一步处理。
主页:
代码语言:javascript复制http://segatalab.cibio.unitn.it/tools/metaphlan/index.html
帮助文档:
代码语言:javascript复制https://github.com/biobakery/MetaPhlAn/wiki/MetaPhlAn-3.0
官方网址:
代码语言:javascript复制http://huttenhower.sph.harvard.edu/metaphlan3
二、软件运行
metaphlan 调用 bowtie2 进行比对,可以输入多种格式结果,例如.fasta, .fastq,bowtie2out 和 sam 等。这里拿之前下载的一个案例 SRS011243 为例。
代码语言:javascript复制cp ./kneaddata_output/SRS011243.denovo_duplicates_marked.trimmed.1_kneaddata_paired_1.fastq SRS011243_1.fq
cp ./kneaddata_output/SRS011243.denovo_duplicates_marked.trimmed.2_kneaddata_paired_1.fastq SRS011243_2.fq
pigz SRS011243_1.fq SRS011243_2.fq
conda activate mpa
#检查下环境中metaphlan_database之前下载的,是否解压缩了
metaphlan --input_type fastq --nproc 12 --bowtie2out metagenome.bowtie2.bz2 SRS011243_1.fq.gz,SRS011243_1.fq.gz -o SRS011243.txt
#还是在下载数据库
Downloading http://cmprod1.cibio.unitn.it/biobakery4/metaphlan_databases/mpa_latest
Downloading http://cmprod1.cibio.unitn.it/biobakery4/metaphlan_databases/mpa_vJan21_CHOCOPhlAnSGB_202103.tar
#中途断了就加上之前下载好的数据库索引如下
--index ~/Software/miniconda3/envs/mpa/lib/python3.7/site-packages/metaphlan/metaphlan_databases/mpa_vJan21_CHOCOPhlAnSGB_202103 #加上刚才下载的数据库索引
#bsub -q fat -n 12 -o %J.log -e %J.err sh meta.sh
三、结果解读
默认输出的结果分为两列,用制表符分割,第一列为物种分类 clades,按照物种很累层级进行排列。每个层级用前缀开头,Kingdom: k__, Phylum: p__, Class: c__, Order: o__, Family: f__, Genus: g__, Species: s__.。第二列是物种对应的相对丰度,不是原始 reads数,所以,每一个层级累积起来总和都为 100%,有些软件输出的是绝对的 reads 数目。
metaphlan 物种分类结果
代码语言:javascript复制#循环处理
#一、首先是kneaddata质控
cat /share/home/xiehs/17.meta/data/hmp/metadata.txt | awk '{if (NR > 1) print $1}' | while read i;do echo "kneaddata -i1 /share/home/xiehs/17.meta/data/hmp/input/${i}/${i}.denovo_duplicates_marked.trimmed.1.fastq -i2 /share/home/xiehs/17.meta/data/hmp/input/${i}/${i}.denovo_duplicates_marked.trimmed.2.fastq -db /share/home/xiehs/17.meta/database/kneadData_databases/human_genome_bowtie2/Homo_sapiens -o ${i}kneaddata_output --remove-intermediate-output -v -t 12 --trimmomatic /share/home/xiehs/Software/miniconda3/envs/biobakery/share/trimmomatic/ --trimmomatic-options 'ILLUMINACLIP:/share/home/xiehs/Software/miniconda3/envs/biobakery/share/trimmomatic/adapters/TruSeq3-PE.fa:2:40:15 SLIDINGWINDOW:4:20 MINLEN:50' --reorder --bowtie2-options '--very-sensitive --dovetail' --run-fastqc-start --run-fastqc-end" >>kne.sh;done;
#echo "parallel -j 4 -a kne.sh" >para.sh
#集群使用32线程
#bsub -q fat -n 48 -o %J.log -e %J.err sh para.sh
#二、跑metaphlan鉴定
cat /share/home/xiehs/17.meta/data/hmp/metadata.txt | awk '{if (NR > 1) print $1}' | while read i;do echo "metaphlan --index ~/Software/miniconda3/envs/mpa/lib/python3.7/site-packages/metaphlan/metaphlan_databases/mpa_vJan21_CHOCOPhlAnSGB_202103 --input_type fastq --nproc 12 /share/home/xiehs/17.meta/illumina/${i}kneaddata_output/${i}.denovo_duplicates_marked.trimmed.1_kneaddata_paired_1.fastq,/share/home/xiehs/17.meta/illumina/${i}kneaddata_output/${i}.denovo_duplicates_marked.trimmed.1_kneaddata_paired_2.fastq --bowtie2out ${i}.bowtie2.bz2 -o ${i}.txt" >>met.sh;done;
#集群使用12线程
#bsub -q fat -n 12 -o %J.log -e %J.err sh met.sh
#结果就是20个txt
#合并
merge_metaphlan_tables.py *txt >../merged_abundance_table.txt
#合并后的物种分类表txt就可以直接去R中处理了。
写在最后:有时间我们会努力更新的。大家互动交流可以前去论坛,地址在下面,复制去浏览器即可访问,弥补下公众号没有留言功能的缺憾。
代码语言:javascript复制bioinfoer.com
有些板块也可以预设为大家日常趣事的分享等,欢迎大家来提建议。