比对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 条评论
登录 后参与评论

相关文章

来自专栏用户2442861的专栏

Python-OpenCV 处理图像(二):滤镜和图像运算

喜欢自拍的人肯定都知道滤镜了,下面代码尝试使用一些简单的滤镜,包括图片的平滑处理、灰度化、二值化等:

511
来自专栏AIUAI

Caffe2 - (二十四) Detectron 之 utils 函数(2)

47511
来自专栏xingoo, 一个梦想做发明家的程序员

推荐系统那点事 —— 基于Spark MLlib的特征选择

在机器学习中,一般都会按照下面几个步骤:特征提取、数据预处理、特征选择、模型训练、检验优化。那么特征的选择就很关键了,一般模型最后效果的好坏往往都是跟特征的选...

2069
来自专栏debugeeker的专栏

《coredump问题原理探究》windows版3.3节函数参数

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/xuzhina/article/detai...

421
来自专栏数据科学与人工智能

【Python环境】Python数据分析入门

本文来分享一下如何通过Python来开始数据分析。 具体内容如下: 数据导入 导入本地的或者web端的CSV文件; 数据变换; 数据统计描述; 假设检验 单样本...

22510
来自专栏算法channel

Spark|有向无环图(DAG)检测

01 — Spark背景介绍 Apache Spark 是专为大规模数据处理而设计的快速通用的计算引擎。Spark 是一种与 Hadoop 相似的开源集群计算环...

3828
来自专栏aCloudDeveloper

LeetCode:152_Maximum Product Subarray | 最大乘积连续子数组 | Medium

题目:Maximum Product Subarray Find the contiguous subarray within an array (contai...

1899
来自专栏wym

python实现opencv学习七:图片色素的数值运算(加减乘除)和逻辑运算(与或非异或)

例图:(若想用下面两张图可另存为图片,若保存的文件无后缀,添加后缀为.jpg即可)

1514
来自专栏开发与安全

VS2008 + Opencv2.1 读取图片像素输出至Excel文件

系统环境: win 7 + VS2008 + Opencv2.1 + Excel 2010 思路:先通过Opencv库函数读取图片存储至IplImage结构体中...

2688
来自专栏北京马哥教育

详解 Python qrcode 二维码模块

1、version:控制二维码的大小,取值范围从1到40。取最小值1时,二维码大小为21*21。取值为 None (默认)或者使用fit=true参数(默认)时...

1000

扫码关注云+社区