老板就是上帝
上帝说要有光,于是就有了光
导师说要有可视化结果,于是我努力把结果可视化~
前情提示
为了知道测序reads中污染物的来源,我们之前学习使用了blast,甚至为了追求速度,我们进一步学习使用了blat和diamond。
做完这些分析之后,我们需要自己去统计,再去可视化(一切都是为了让老板的体验更佳)。
但是,一次寻找public data的过程中,我偶然间发现NCBI对于上传的数据也有自己的taxonomy analysis。NCBI使用一个名为NCBI SRA Taxonomy Analysis Tool (STAT)的官方软件进行分析,其结果也可以直接用于Krona进行可视化
接下来,我们随便找个数据分别看一下软件的运行结果以及可视化结果
其中可视化结果中,从内到外,分别是界门纲目科属种。而绿色部分表示那些没有识别出的部分。
STAT (NCBI SRA Taxonomy Analysis Tool)
对于STAT的运行结果,它给出了数据中reads在不同分类节点上的百分比,若一个reads可以同时比对到超过两个分类节点,则这个reads被划分到分类层级最低的共享分类中。比如,当一个read映射到属于同一属的两个物种时,它被分配在属这个分类层级。
当一个reads可以比对到一个物种的多个子分类时,那么高分类层级的占比通常高于低层级分类的占比。
STAT的分类方法,类似于blast,将参考基因组库,按照k-mer建立索引,然后将reads进行比对/query。接着将比对上的reads进行分类,计算百分比。
STAT的下载安装相当复杂(似乎尚不支持conda安装),以下是具体代码:
mkdir ~/SRC ; cd ~/SRC
git clone https://github.com/ncbi/ngs.git
git clone https://github.com/ncbi/ncbi-vdb.git
git clone https://github.com/ncbi/ngs-tools.git
cd ngs/ngs-sdk ; ./configure -p=~/NCBI && make
cd ../../ncbi-vdb ; ./configure -p=~/NCBI && make
cd ../ngs-tools ; ./configure -p=~/NCBI
具体用法可以参考示例文件,在 examples
文件夹里有一些shell脚本
参考资料
Krona
Krona可以让我们对分类数据进行可视化,得到的数据结果可以放大缩小,既可以通过excel插件绘制可视化结果也可以使用KronaTools软件包。它生成的可视化图是交互式的,可以自行调整展示方式。
下载安装
这里我们以KronaTools为例介绍Krona的使用,首先我们下载软件包:
wget -c https://github.com/marbl/Krona/archive/v2.7.1.tar.gz
tar -zxvf v2.7.1.tar.gz
cd /Krona-2.7.1/KronaTools
mkdir taxonomy
./install.pl --prefix $PWD --taxonomy $PWD
# 此处在--prefix后不需要加bin,会自动在路径下生成一个bin目录
-
--prefix
指定bin目录所在的路径 -
--taxonomy
指定updateTaxonomy.sh
或updateAccessions.sh
这两个脚本运行时,得到的taxonomy(物种分类)文件所存储的位置
Krona中的一些工具依赖于NCBI taxonomy,要想正确使用,需要运行当前路径下的 updateTaxonomy.sh
脚本。该脚本下载一个本地的taxonomy数据库,不到100M。也可以使用该脚本进行NCBI taxonomy数据库的更新。
此外,还有一些工具依赖于accession 和 taxonomy ID之间的转换查找,因此需要运行 updateAccession.sh
脚本,安装本地的accession-to-taxonomyID数据库,该数据库较大,约16GB。
除了在安装过程中,进行数据库的下载外,也可以进行手动下载,后续再进行配置。具体方式见Krona的github(https://github.com/marbl/Krona/wiki/Installing)
具体使用
Krona支持多种输入文件格式,我最常用blast进行比对,因此以blast的输出结果进行举例。在Krona的安装目录下存在一个 ktImportBLAST
脚本,它可以直接利用blast输出结果产生可视化网页。具体脚本如下
cd ./bin
./ktImportBLAST /test_path/test.blast.out -o /test_path/test.blast.html
-
-o
参数指定输出的结果文件名称,若需要分类的reads数目过多,则Krona会新产生一个同名文件夹存放这些额外信息。
其他的输入文件格式此处不再展开讲,具体可参考软件github主页(https://github.com/marbl/Krona/wiki/KronaTools)