BLAST结果可视化

老板就是上帝

上帝说要有光,于是就有了光

导师说要有可视化结果,于是我努力把结果可视化~

前情提示

为了知道测序reads中污染物的来源,我们之前学习使用了blast,甚至为了追求速度,我们进一步学习使用了blat和diamond。

做完这些分析之后,我们需要自己去统计,再去可视化(一切都是为了让老板的体验更佳)。

但是,一次寻找public data的过程中,我偶然间发现NCBI对于上传的数据也有自己的taxonomy analysis。NCBI使用一个名为NCBI SRA Taxonomy Analysis Tool (STAT)的官方软件进行分析,其结果也可以直接用于Krona进行可视化

接下来,我们随便找个数据分别看一下软件的运行结果以及可视化结果

Taxonomy analysis

Krona visualize

其中可视化结果中,从内到外,分别是界门纲目科属种。而绿色部分表示那些没有识别出的部分。

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.shupdateAccessions.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

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 199,340评论 5 467
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 83,762评论 2 376
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 146,329评论 0 329
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,678评论 1 270
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,583评论 5 359
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 47,995评论 1 275
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,493评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,145评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,293评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,250评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,267评论 1 328
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,973评论 3 316
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,556评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,648评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,873评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,257评论 2 345
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 41,809评论 2 339

推荐阅读更多精彩内容