三代测序宏基因组物种分类鉴定

一、mock 数据集

人为添加 10 种微生物，其中包括 8 株细菌，两株真菌。分为两种模式，一种按比例平均分配，称为 Even 数据集，8 株细菌各占 8%，2 株真菌各占 4%。另一种按照对数进行分配，称为 Log 数据集。并且包括 illumina 与 nanopore 数据。

https://github.com/LomanLab/mockcommunity
https://lomanlab.github.io/mockcommunity/

二、土壤样品

文章列表：https://www.biorxiv.org/content/10.1101/2020.04.08.032540v1
练习数据：https://www.ebi.ac.uk/ena/browser/view/PRJEB36155

2012 年加拿大 Northern Alberta 油砂尾矿池附近海藻细菌培养分离样品，2019 年重新培养提取。

脚本路径：

https://zenodo.org/record/3745531#.Y1lAd7ZBxPa

《Complete and validated genomes from a metagenome》

数据一般都在文章结尾的“Data availability”部分，从中找到 BioProject 号或者 SRA 号即可。例如该文章中给出了数据的 BioProject 号为 PRJEB36115。

三、centrifuge 物种鉴定

centrifuge 的使用非常简单，输入数据包含测序的数据以及索引文件。可支持二代和三代测序数据，输入为 fastq 格式文件即可，也支持 fasta 格式以及原始 qseq 格式文件，同时支持pairend 数据，也支持压缩格式。其中索引只写前缀名即可。

#centrifuge 进行物种分类鉴定
centrifuge -x centrifuge_h+p+v_20200318/hpv -U nanopore.fastq.gz --report-file
report.tsv -S result.tsv -p 64 >centrifuge.log

四、结果解读

centrifuge 默认会输出两个文件，分别是按照 reads 进行统计的结果与按照物种进行统计的结果。

1、按照 reads 进行统计的结果 centrifuge_output.tsv

centrifuge 结果展示

该文件一共分类 8 列 。

1：原始 read ID ;

2、比对到数据库中的序列 ID，如果使用的是 Refseq 数据库或者 nt 库，则是序列的 AccessionID;

3、物种分类 ID，第二列比对上序列对应的物种分类 ID;

4、classification 的分值，比对上的序列之和;

5、第二好比对结果分值;

6、比对到序列的长度；

7、比对的 reads 长度；

8、这条 reads 比对上多少个物种序列；

2、按照比对上的物种进行统计 centrifuge_report.tsv

1、比对上物种名字，如果鉴定不到种，则上升一级；

2、物种分类 ID；

3、物种分类层级 rank；

4、对应基因组大小；

5、比对到的 reads 数目，包括多重比对的结果；

6、唯一比对上的 reads 数目；

7、比对的丰度，比对上区域/基因组长度。

五、过滤结果

由于序列相似性的缘故，一条序列可能会比对到数据库中多个物种，Centrifuge 原始的结果会鉴定到很多物种，这就需要对原始数据进行过滤，通常选择每条序列最优的比对。然后根据每个物种比对上的 reads 数进行过滤，同时也可以根据鉴定到的物种水平进行筛选。

awk -F "\t" '{if ($3=="species" && $6 >5) print $1"\t"$6}' 0.01_report.tsv >0.01.txt

当然也可以用R语言去筛选和排序表格。

写在最后：有时间我们会努力更新的。大家互动交流可以前去论坛，地址在下面，复制去浏览器即可访问，弥补下公众号没有留言功能的缺憾。

bioinfoer.com

有些板块也可以预设为大家日常趣事的分享等，欢迎大家来提建议。