比对NR库看看物种分布【直播】我的基因组88

前面我提前了我的基因组测序数据里面的未成功比对到人类基因组上面的那些fastq序列,也用了软件把它们组装成fasta序列,这些序列的功能是未知的,可以通过比对到NCBI的NT/NR库来给他们注释一下。

NR库是Non-redundant protein sequences from GenPept, Swissprot, PIR, PDF, PDB, and NCBI RefSeq,得去ftp://ftp.ncbi.nih.gov/blast/db/ 下载所有gz结尾的文件,并且解压到同一个目录即可。

最终还是得做这件事了,之前懒得做就是因为NCBI为BLAST提供的NR库实在是太大了。但是我无意中看到同一个服务器的朋友就在跑nr库的blast,就直接借用他的啦,可以看到解压后有159G,包括70G的fasta序列,如下:

$cd /home/chenhuang/database/NCBI/NR$ du -h159G    .$ ls -lh |tail |cut -d" " -f 5-10 22M May 31 17:18 nr.46.pin957M May 31 17:18 nr.46.psq502M May 31 17:20 nr.47.phr 21M May 31 17:20 nr.47.pin957M May 31 17:20 nr.47.psq 83M May 31 17:20 nr.48.phr4.7M May 31 17:20 nr.48.pin228M May 31 17:20 nr.48.psq

虽然不需要自己下载nr库,也不需要自己建库,但是blast软件还是要下载的。我以前写过说明书:NCBI的blast++软件使用说明书 当然,如果你没有朋友帮你下载,你还是得自己来。

cd ~/biosoft mkdir blast && cd blast wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/LATEST/ncbi-blast-2.6.0+-x64-linux.tar.gztar ncbi-blast-2.6.0+-x64-linux.tar.gz 

比对分为好几种,blastn, blastp,blastx,tblastn,运行的方法都是一样的

blastp -query seq.fasta -out seq.blast -db dbname -outfmt 6 -evalue 1e-5 -num_descriptions 10 -num_threads 8 

blastp:

protein query

->

protein sequence database

blastn:

nucleotide query

->

nucleotide sequence database

tblastn:

protein query

->

nucleotide sequence database

blastx:

nucleotide query

->

protein sequence database

blastp:将待查询的蛋白质序列及其互补序列一起对蛋白质序列数据库进行查询;blastn:将待查询的核酸序列及其互补序列一起对核酸序列数据库进行查询;blastx:先将待查询的核酸序列按六种可读框架(逐个向前三个碱基和逐个向后三个碱基读码)翻译成蛋白质序列,然后将翻译结果对蛋白质序列数据库进行查询;tblastn:先将核酸序列数据库中的核酸序列按六种可读框架翻译成蛋白质序列,然后将待查询的蛋白质序列及其互补序列对其翻译结果进行查询;tblastx:先将待查询的核酸序列和核酸序列数据库中的核酸序列按六种可读框架翻译成蛋白质序列,然后再将两种翻译结果从蛋白质水平进行查询。

参数说明:

  • -query: 输入文件路径及文件名
  • -out:输出文件路径及文件名
  • -db:格式化了的数据库路径及数据库名
  • -outfmt:输出文件格式,总共有12种格式,6是tabular格式对应BLAST的m8格式
  • -evalue:设置输出结果的e-value值
  • -num_descriptions:tabular格式输出结果的条数
  • -num_threads:线程数

那么我的代码是:

cd ~/data/project/myGenome/gatk/jmzeng/unmapped/~/biosoft/blast/ncbi-blast-2.6.0+/bin/blastn -query output_prefix.contigs.fa \-out seq.blast -db /home/chenhuang/database/NCBI/NR/nr \-outfmt 6 -evalue 1e-5  -num_threads 5

输出结果解释

重点是-outfmt 6,也就是之前版本的m 8格式

结果中从左到右每一列的意义分别是:

  • [00] Query id
  • [01] Subject id
  • [02] % identity
  • [03] alignment length
  • [04] mismatches
  • [05] gap openings
  • [06] q. start
  • [07] q. end
  • [08] s. start
  • [09] s. end
  • [10] e-value
  • [11] bit score

简单看前几行信息,如下:

2__len__515    XP_015633936.1  100.000 127 0   0   510 130 155 281 1.13e-59    1962__len__515    CAH68354.1  99.213  127 1   0   510 130 155 281 5.18e-59    1942__len__515    EEC78288.1  99.213  127 1   0   510 130 206 332 2.88e-58    1942__len__515    XP_015691624.1  79.688  128 13  3   510 130 95  210 1.01e-44    1552__len__515    BAH92880.1  100.000 85  0   0   283 29  209 293 3.52e-25    1082__len__515    CAE03461.1  100.000 85  0   0   283 29  433 517 1.40e-24    1082__len__515    OEL26449.1  72.277  101 18  3   513 217 76  168 2.68e-23    1002__len__515    XP_003579498.1  69.072  97  23  3   507 229 177 270 3.79e-23    1022__len__515    BAJ88079.1  63.636  99  28  4   507 217 158 250 1.38e-16    84.72__len__515    XP_020196607.1  64.646  99  27  4   513 223 157 249 2.10e-16    84.3

可以看到第二列的NR库里面的蛋白ID仍然是多种多样的,虽然是以refseq ID为主,需要把它们对应到各个物种,这样才知道我们的序列的物种分布。而且这个比对到NR库也实在是太耗费时间了,整整24个小时才处理了667条序列。

blast结果的详细比对结果。注意比对到的序列长度。评价一个blast结果的标准主要有三项,E值(Expect),一致性(Identities),缺失或插入(Gaps)。加上长度的话,就有四个标准了。

E值(Expect):表示随机匹配的可能性,例如,E=1,表示在目前大小的数据库中,完全由机会搜到对象数的平均值为1.E值越大,随机匹配的可能性也越大。E值接近零或为零时,具本上就是完全匹配了。通常来讲,我们认为E值小于10-5 就是比较可性的S值结果。我们可以想象,相同的数据库,E=0.001时如果有1000条都有机会S值比现在这个要高的话,那么不E设置为10-6时可能就会只得到一条结果,就是S值最可靠的那个。但是E值也不是万能的。它在以下几个情况下有局限性:

1)当目标序列过小时,E值会偏大,因为无法得到较高的S值。2)当两序列同源性虽然高,但有较大的gap(空隙)时,S值会下降。这个时候gap scores就非常有用。3)有些序列的非功能区有较低的随机性时,可能会造成两序列较高的同源性。

E值总结:

E值适合于有一定长度,而且复杂度不能太低的序列。当E值小于10-5时,表明两序列有较高的同源性,而不是因为计算错误。当E值小于10-6时,表时两序列的同源性非常高,几乎没有必要再做确认。

一致性(Identities):或相似性。匹配上的碱基数占总序列长的百分数。

Score得分值越高说明同源性越好;Expect期望值越小比对结果越好,说明因某些原因而引起的误差越小;Identities是同源性(相似性),例中所示比对的1299个碱基中只有35个不配,其他97%相同;

Gaps是指多出或少的碱基或缺失的碱基数;缺失或插入(Gaps):插入或缺失。用"—"来表示。

Strand=plus/plus指两条序列方向相同,如果是plus/minus,即意味着一条是5'到3',一条是3'到5',或一条是正向,另一条是反向序列。

这些ID的具体描述信息,都是那个70G的FASTA序列数据里面。

grep '^>' nr |cut -d"]" -f 1  |head>WP_003131952.1 30S ribosomal protein S18 [Lactococcus lactis>XP_642131.1 hypothetical protein DDB_G0277827 [Dictyostelium discoideum AX4>XP_642837.1 hypothetical protein DDB_G0276911 [Dictyostelium discoideum AX4>WP_000184067.1 MULTISPECIES: MbtH family protein [Bacillus>WP_007051162.1 MULTISPECIES: argininosuccinate lyase [Bifidobacterium>WP_000135199.1 MULTISPECIES: 30S ribosomal protein S18 [Bacteria>WP_003251213.1 MULTISPECIES: leucyl/phenylalanyl-tRNA--protein transferase [Pseudomonas>WP_003409891.1 MULTISPECIES: SecB-like chaperone [Mycobacterium tuberculosis complex>WP_000379821.1 MULTISPECIES: O-acetyltransferase OatA [Staphylococcus>WP_000332037.1 MULTISPECIES: ribonucleoside-diphosphate reductase 1 subunit beta [Proteobacteria

这里我只关心物种,就不需要知道这个蛋白的意义了,简单的对应一下。简单的脚本制作ID对应文件,4.4G,1.23亿条序列。不过结果比较诡异

      2 Anas platyrhynchos      2 Meleagris gallopavo      2 Zea mays      3 Homo sapiens      3 Nasonia vitripennis      3 Oryza brachyantha      9 Gallus gallus     10 Oryza sativa Indica Group     11 Arabidopsis thaliana     77 Oryza sativa Japonica Group

可能是我想的太简单了,大部分是水稻,拟南芥的序列,还有玉米什么的,总之都是植物。

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-09-12

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏MyBlog

软件测试方法课程笔记(2)

为了产生少量的测试用例, 并且可以测试大部分的情况, 我们可以使用等价类划分的方法 比如对于输入值是范围值, 我们可以使用等价类划分成范围内的和不是范围内的两...

1242
来自专栏Python小屋

猜数游戏用Python应该这样写

from random import randint def guess(): #随机生成一个整数 value = randint(1,1000...

2789
来自专栏钱坤的专栏

cache 淘汰算法:LIRS 算法

LIRS 算法是非常优秀的 cache 淘汰算法,被用于 mysql 5.1之后的版本,这篇文章主要来源于对 LIRS 的发表论文的翻译。

1.6K3
来自专栏企鹅号快讯

从图灵机开始

说到图灵机,我们首先要说说图灵这个人。笔者觉得我们这种搞计算机的人都应该知道并记得这个人。 图灵,1912年6月23日生于英国帕丁顿。是数学家、密码破译专家,当...

2068
来自专栏小樱的经验随笔

从入门到精通之Boyer-Moore字符串搜索算法详解

本文讲述的是Boyer-Moore算法,Boyer-Moore算法作为字符串搜索算法,兴趣之下就想了解这个算法,发现这个算法一开始还挺难理解的,也许是我理解能力...

3428
来自专栏生信技能树

【直播】我的基因组73:在IGV看看indel是啥样子

前面我们特意用scalpel软件来找indel,期待它会有一些出彩的表现,当然我还没来得及比较它找到的INDEL跟GATK等工具区别在哪里,不过我们先在IGV里...

4779
来自专栏pydata

Matlab C混合编程

在MATLAB中可调用的C或Fortran语言程序称为MEX文件。MATLAB可以直接把MEX文件视为它的内建函数进行调用。MEX文件是动态链接的子例程,MAT...

912
来自专栏程序员互动联盟

二维码是如何实现的?

二维条码是指在一维条码的基础上扩展出另一维具有可读性的条码,使用黑白矩形图案表示二进制数据,被设备扫描后可获取其中所包含的信息。一维条码的宽度记载着数据,而其长...

3255
来自专栏腾讯AlloyTeam的专栏

教你用 webgl 快速创建一个小世界

Webgl的魅力在于可以创造一个自己的3D世界,但相比较canvas2D来说,除了物体的移动旋转变换完全依赖矩阵增加了复杂度,就连生成一个物体都变得很复杂……这...

1.7K0
来自专栏生信宝典

R包ggseqlogo 绘制seq logo图

在生物信息分析中,经常会做序列分析图(sequence logo),这里的序列指的是核苷酸(DNA/RNA链中)或氨基酸(在蛋白质序列中)。sequence l...

1943

扫码关注云+社区