宏基因组数据分析–如何完全去除人源序列–计算微生物丰度

下机得到原始fastq文件。经过fastp软件处理得到质量相对较高的cleandata,之后再使用kraken2进行物种注释。

1.如何完全去除人源序列

上一篇文章中讲到去除人源序列通过bowtie2(蝴蝶结)比对到hg38基因组，利用samtools去除比对上的基因，但是去除的人源总是不干净，还是含有40%的人基因，（因为人基因组上的SNP和snv是普遍存在的，所以会有残留）

如果不将人的去除干净，你再去计算丰度（把人的基因也算进去了）是没有意义的

今天捣鼓了半天，终于找到了好的方法：

利用krakentools里的 extract_kraken_reads.py 

先进行kraken2，得到kraken2比对的详细文件 xx.kraken(不是report)

再利用这个脚本：

extract_kraken_reads.py -k xxx.kraken -s1 read1.fq -s2 reads2.fq -o extracted1.fq -o2 extracted2.fq --taxid 9606  --exclude

taxid 9606 是在NCBI中智慧人类的taxonomy IDs

2.计算丰度

下面是一篇2020年的一篇文章

Kraken是用于分类宏基因组测序数据的非常快速，准确的程序。它需要一组读数，重叠群或其他DNA序列，并为每一个分配一个分类标签（物种，属等）。发现有些直接使用Kraken进行丰量估计 - 估计样本中物种的相对比例 - 并基于假设Kraken的输出可以用于这种方式发表论文。但是，这是不正确的。如果您给Kraken一套宏基因组阅读器进行分类，它将为每个阅读分配最具体的标签。但是很多时候，这些标签并不属于物种层面。例如，如果150bp的读数与两个不同的物种完全相同，Kraken将把它分配给它们的最低共同祖先（LCA），它可能属于等级或更高。对于含有两个或更多高度相似物种的样品，这意味着物种特异性读数的数量可能远远少于预期。（我们应该注意到，Kraken经常在应变级别指定读取。）

为了解决这个问题，我们开发了Bracken：用KrakEN分类后的贝叶斯重估丰度。Bracken使用贝叶斯算法和Kraken分类结果来估计宏基因组样本的物种水平或属水平丰度。这个新工具在我们的PeerJ论文或Bracken软件网站中有详细描述。

Bracken将Kraken的输出转化为物种丰度的估计。这在同一样品中出现高度相似的物种时特别有用。例如，牛分枝杆菌和结核分枝杆菌的基因组是99.95％相同的。因此，Kraken将这些物种的绝大多数分类为分枝杆菌属（因为它们是无法区分的）。但是，特定样本中的一些读段通常来自基因组的独特（或物种特异性）部分。Bracken使用这些信息，加上关于姊妹物种之间相似性的信息，将属性级别的读数“下推”到物种级别。如果您有兴趣了解样本中每种分枝杆菌种类的存在情况，可以使用Kraken，然后使用Bracken来估算每种物种的相对丰度。

那么现在能直接拿kraken2的分析数据直接做丰度分析吗？

                                这是2020年kraken的版本

在后面的版本更新中并没有发现对此问题进行升级

下面拿这两个软件的报告进行对比一下

所以最终看丰度报告的话，还是建议bracken!!!!!!

bracken -d KRAKEN_DB -i xxx.kreport -o xxx.bracken -r 300 -l S 

-r是序列长度（根据测序时确定，单端还是双端）

-l S(确定分类水平，具体自行看bracken说明

得到的文件xxx.bracken是可以用Excel打开的，大家可以发现已经完全没有人源的了，完全反映了微生物的丰度情况，可以操作Excel自行按照丰度进行排序(现在丰度较高的微生物已经能达到20左右了)