比对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可能是我想的太简单了，大部分是水稻，拟南芥的序列，还有玉米什么的，总之都是植物。