山内のサイト

Blastのダウンロード †

Blastのダウンロード　Download BLAST Software and Databases Documentation

wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.9.0+-x64-linux.tar.gz
tar -xvf ncbi-blast-2.9.0+-x64-linux.tar.gz

インストール　Installation - BLAST; Command Line Applications User Manual - NCBI Bookshelf

ここでは　/usr/local/bin/ncbi-blast-2.9.0+/　に展開されているので、

ln -s /usr/local/ncbi-blast-2.9.0+ blast

で /usr/local/blast を作っておき、PATHに /usr/local/blast/bin を追加。

使い方のマニュアルは　BLAST® Command Line Applications User Manual - NCBI Bookshelf にある。

BLAST configuration fileを作るのが良いらしい。

↑

自分用のデータベースを作る †

FASTAから作れるらしい。Building a BLAST database with local sequences - BLAST® Command Line Applications User Manual - NCBI Bookshelf サンプルで試してみる。

>seq1
MSFSTKPLDMATWPDFAALVERHNGVWGGCWCMAFHAKGSGAVGNREAKEARVREGSTHAALVFDGSACVGWCQFGPTGE
LPRIKHLRAYEDGQAVLPDWRITCFFSDKAFRGKGVAAAALAGALAEIGRLGGGTVESYPEDAQGRTVAGAFLHNGTLAM
>seq2
MKAIDLKAEEKKRLIEGIQDFFYEERNEEIGIIAAEKALDFFLSGVGKLIYNKALDESKIWFSRRLEDISLDYELLYK
>seq3 
MTLAAAAQSATWTFIDGDWYEGNVAILGPRSHAMWLGTSVFDGARWFEGVAPDLELHAARVNASAIALGLAPNMTPEQIV
GLTWDGLKKFDGKTAVYIRPMYWAEHGGYMGVPADPASTRFCLCLYESPMISPTGFSVTVSPFRRPTIETMPTNAKAGCL
YPNNGRAILEAKARGFDNALVLDMLGNVAETGSSNIFLVKDGHVLTPAPNGTFLSGITRSRTMTLLGDYGFRTTEKTLSV
RDFLEADEIFSTGNHSKVVPITRIEGRDLQPGPVAKKARELYWDWAHSASVG
>seq4
MRSFFHHVAAADPASFGVAQRVLTIPIKRAHIEVTHHLTKAEVDALIAAPNPRTSRGRRDRTFLLFLARTGARVSEATGV
NANDLQLERSHPQVLLRGKGRRDRVIPIPQDLARALTALLAEHGIANHEPRPIFIGARQERLTRFGATHIVRRAAAQAVT
IKPALAHKPISPHIFRHSLAMKLLQSGVDLLTIQAWLGHAQVATTHRYAAADVEMMRKGLEKAGVSGDLGLRFRPNDAVL
QLLTSI
>seq5
MTISRVCGSRTEAMLTNGQEIAMTSILKSTGAVALLLLYTLTANATSLMISPSSIERVAPDRAAVFHLRNQMDRPISIKV
RVFRWSQKGGVEKLEPTGDVVASPISAQLSPNGNRAVRVVRVSKEPLRSEEGYRVVIDEADPTRNTPEAESLSARHVLPV
LFRPPDVLGPEIELSLTRSDGWLMLVVENKGASRLRRSDVTLAQGSAGIARREGFVGYVLPGLTRHWRVGREDSYSGGIV
TVSANSSGGAIGEQLVVSGR
>seq6
TTLLLQVPIGWGVLHQGGALVVLGFAIAHWRGFVGTYTRDTAIEMRD

これに対して

makeblastdb -in test.fsa -parse_seqids -blastdb_version 5 -title "Cookbook demo" -out test -dbtype prot

Building a new DB, current time: 09/03/2019 14:35:54
New DB name:   /home/yamanouc/src/Kishimoto2017/test-Blast/test.fsa
New DB title:  Cookbook demo
Sequence type: Protein
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 6 sequences in 0.0214489 seconds.

更に、比較入力用の配列として test.fasta を予め（データベースデータに使ったseq1だけを含むファイルとして）作っておく。

これを使って、

blastp -db test < test.fasta

とすると、動作する。

BLASTP 2.9.0+


Reference: Stephen F. Altschul, Thomas L. Madden, Alejandro A.
Schaffer, Jinghui Zhang, Zheng Zhang, Webb Miller, and David J.
Lipman (1997), "Gapped BLAST and PSI-BLAST: a new generation of
protein database search programs", Nucleic Acids Res. 25:3389-3402.


Reference for composition-based statistics: Alejandro A. Schaffer,
L. Aravind, Thomas L. Madden, Sergei Shavirin, John L. Spouge, Yuri
I. Wolf, Eugene V. Koonin, and Stephen F. Altschul (2001),
"Improving the accuracy of PSI-BLAST protein database searches with
composition-based statistics and other refinements", Nucleic Acids
Res. 29:2994-3005.



Database: Cookbook demo
           6 sequences; 1,083 total letters



Query= seq1

Length=160
                                                                     Score        E
Sequences producing significant alignments:                          (Bits)     Value

seq1 unnamed protein product                                          328        2e-121
seq4 unnamed protein product                                          14.2       6.5   
seq3 unnamed protein product                                          14.2       6.9   


>seq1 unnamed protein product
Length=160

 Score = 328 bits (840),  Expect = 2e-121, Method: Compositional matrix adjust.
 Identities = 160/160 (100%), Positives = 160/160 (100%), Gaps = 0/160 (0%)

Query  1    MSFSTKPLDMATWPDFAALVERHNGVWGGCWCMAFHAKGSGAVGNREAKEARVREGSTHA  60
            MSFSTKPLDMATWPDFAALVERHNGVWGGCWCMAFHAKGSGAVGNREAKEARVREGSTHA
Sbjct  1    MSFSTKPLDMATWPDFAALVERHNGVWGGCWCMAFHAKGSGAVGNREAKEARVREGSTHA  60

Query  61   ALVFDGSACVGWCQFGPTGELPRIKHLRAYEDGQAVLPDWRITCFFSDKAFRGKGVAAAA  120
            ALVFDGSACVGWCQFGPTGELPRIKHLRAYEDGQAVLPDWRITCFFSDKAFRGKGVAAAA
Sbjct  61   ALVFDGSACVGWCQFGPTGELPRIKHLRAYEDGQAVLPDWRITCFFSDKAFRGKGVAAAA  120

Query  121  LAGALAEIGRLGGGTVESYPEDAQGRTVAGAFLHNGTLAM  160
            LAGALAEIGRLGGGTVESYPEDAQGRTVAGAFLHNGTLAM
Sbjct  121  LAGALAEIGRLGGGTVESYPEDAQGRTVAGAFLHNGTLAM  160


>seq4 unnamed protein product
Length=246

 Score = 14.2 bits (25),  Expect = 6.5, Method: Compositional matrix adjust.
 Identities = 8/27 (30%), Positives = 11/27 (41%), Gaps = 0/27 (0%)

Query  101  RITCFFSDKAFRGKGVAAAALAGALAE  127
            R+T F +    R     A  +  ALA 
Sbjct  141  RLTRFGATHIVRRAAAQAVTIKPALAH  167


>seq3 unnamed protein product
Length=292

 Score = 14.2 bits (25),  Expect = 6.9, Method: Compositional matrix adjust.
 Identities = 6/11 (55%), Positives = 8/11 (73%), Gaps = 0/11 (0%)

Query  56   GSTHAALVFDG  66
            GS++  LV DG
Sbjct  192  GSSNIFLVKDG  202



Lambda      K        H        a         alpha
   0.321    0.136    0.447    0.792     4.96 

Gapped
Lambda      K        H        a         alpha    sigma
   0.267   0.0410    0.140     1.90     42.6     43.6 

Effective search space used: 117390


  Database: Cookbook demo
    Posted date:  Sep 3, 2019  2:57 PM
  Number of letters in database: 1,083
  Number of sequences in database:  6



Matrix: BLOSUM62
Gap Penalties: Existence: 11, Extension: 1
Neighboring words threshold: 11
Window for multiple hits: 40

出力を細かく制御しようとすると、たとえば

blastp -db test -outfmt "10 qseqid sseqid pident length mismatch evalue bitscore" -out test.output < test.fasta

とすると、

seq1,seq1,100.000,160,0,2.25e-121,328
seq1,seq4,29.630,27,19,6.5,14.2
seq1,seq3,54.545,11,5,6.9,14.2

が得られる。詳細はblastp -helpおよびDisplay BLAST search results with custom output format - BLAST® Command Line Applications User Manual - NCBI Bookshelfを参照。

Evalueの値の制限を緩めて（大きくして）、もっと幅広く一致を取ると、（ここでは入力を元のtest.fsaの全配列にしている）

blastp -db test -outfmt "10 qseqid sseqid pident length mismatch evalue bitscore" -evalue 100 < test.fsa

だと、ほぼ全部の組合せについて出力される。

seq1,seq1,100.000,160,0,2.25e-121,328
seq1,seq4,29.630,27,19,6.5,14.2
seq1,seq3,54.545,11,5,6.9,14.2
seq1,seq3,25.000,40,30,22,12.7
seq1,seq3,50.000,8,4,28,12.3
seq1,seq3,33.333,21,14,40,11.9
seq1,seq3,33.333,9,6,57,11.5
seq1,seq5,25.000,20,10,44,11.9
seq1,seq5,35.000,20,12,90,10.8
seq2,seq2,100.000,78,0,5.31e-54,151
seq2,seq3,33.333,12,8,56,10.4
seq3,seq3,100.000,292,0,0.0,602
seq3,seq5,19.403,67,39,0.28,19.6
seq3,seq5,28.000,25,18,2.5,16.5
seq3,seq5,29.730,37,25,4.9,15.8
seq3,seq5,50.000,8,4,59,12.3
seq3,seq4,27.711,83,39,8.1,15.0
seq3,seq4,66.667,6,2,17,13.9
seq3,seq4,27.368,95,49,45,12.7
seq3,seq1,54.545,11,5,14,14.2
seq3,seq1,25.000,40,30,33,13.1
seq3,seq1,50.000,8,4,50,12.3
seq3,seq1,33.333,9,6,99,11.5
seq4,seq4,100.000,246,0,0.0,489
seq4,seq3,28.421,95,48,2.2,16.5
seq4,seq3,27.473,91,45,6.0,15.0
seq4,seq3,66.667,6,2,13,14.2
seq4,seq5,38.462,13,8,8.9,14.6
seq4,seq5,42.857,21,9,35,12.7
seq5,seq5,100.000,260,0,0.0,519
seq5,seq3,19.403,67,39,0.26,19.6
seq5,seq3,28.000,25,18,2.6,16.5
seq5,seq3,29.730,37,25,4.7,15.4
seq5,seq3,50.000,8,4,51,12.3
seq5,seq4,38.462,13,8,10.0,14.6
seq5,seq4,42.857,21,9,32,13.1
seq5,seq1,25.000,20,10,77,11.5
seq6,seq6,100.000,47,0,6.34e-33,95.9
↑

本番データでやってみる †

実際のデータでやってみる。AE000512側のCDSのfastaファイルをAE000512CDS.fastaとして作ってあるので、これをデータベースに変換する。

makeblastdb -in AE000512CDS.fasta -parse_seqids -blastdb_version 5 -title "Thermo" -out AE000512 -dbtype nucl


Building a new DB, current time: 09/03/2019 15:23:58
New DB name:   /home/yamanouc/src/Kishimoto2017/test-Blast/AE000512
New DB title:  Thermo
Sequence type: Nucleotide
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 1846 sequences in 0.117031 seconds.

このデータベースを使って、blastn を起動し、AP012030側のCDSのfastaファイルAP012030CDS.fastaを入力として与える。

blastn -db AE000512 -outfmt "10 qseqid sseqid pident length mismatch evalue bitscore" -out AP012030blast.out -evalue 100 < AP012030CDS.fasta

としたところ、得られたのは

EcoliECDH1ME8569_1662,ThermoTM_0740,100.000,30,0,3.16e-08,56.5

だけであった。Evalueをここまで広げても１つだけだった。

逆を試してみよう。つまり、データベースにAP012030を置き、blast入力にAE000512を 入れてみる。

データベースを作る。

makeblastdb -in AP012030CDS.fasta -parse_seqids -blastdb_version 5 -title "Ecoli" -out AP012030 -dbtype nucl

AE000512CDS.fastaと比較する。

blastn -db AP012030 -outfmt "10 qseqid sseqid pident length mismatch evalue bitscore" -out AE000512blast.out -evalue 100 < AE000512CDS.fasta

結果は同じで、

ThermoTM_0740,EcoliECDH1ME8569_1662,100.000,30,0,7.23e-08,56.5

であった。

なお、この２つの遺伝子は、AP012030側が

    gene            complement(1784023..1785951)
                    /gene="thrS"
                    /locus_tag="ECDH1ME8569_1662"
                    /gene_synonym="b1719"
                    /gene_synonym="ECK1717"
                    /gene_synonym="JW1709"
    CDS             complement(1784023..1785951)
                    /gene="thrS"
                    /locus_tag="ECDH1ME8569_1662"
                    /gene_synonym="b1719"
                    /gene_synonym="ECK1717"
                    /gene_synonym="JW1709"
                    /function="Threonyl-tRNA synthetase"
                    /codon_start=1
                    /transl_table=11
                    /product="threonyl-tRNA synthetase"
                    /protein_id="BAJ43518.1"

配列が

ATGCCTGTTATAACTCTTCCTGATGGCAGCCAACGCCATTACGATCACGCTGTAAGCCCCATGGATGTTGCGCTGGACATTGGTCCAGGTCTGGCGAAAGCCTGTATCGCAGGGCGCGTTAATGGCGAACTGGTTGATGCTTGCGATCTGATTGAAAACGACGCACAACTGTCGATCATTACCGCCAAAGACGAAGAAGGTCTGGAGATCATTCGTCACTCCTGTGCGCACCTGTTAGGGCACGCGATTAAACAACTTTGGCCGCATACCAAAATGGCAATCGGCCCGGTTATTGACAACGGTTTTTATTACGACGTTGATCTTGACCGCACGTTAACCCAGGAAGATGTCGAAGCACTCGAGAAGCGGATGCATGAGCTTGCTGAGAAAAACTACGACGTCATTAAGAAGAAAGTCAGCTGGCACGAAGCGCGTGAAACTTTCGCCAACCGTGGGGAGAGCTACAAAGTCTCCATTCTTGACGAAAACATCGCCCATGATGACAAGCCAGGTCTGTACTTCCATGAAGAATATGTCGATATGTGCCGCGGTCCGCACGTACCGAACATGCGTTTCTGCCATCATTTCAAACTAATGAAAACGGCAGGGGCTTACTGGCGTGGCGACAGCAACAACAAAATGTTGCAACGTATTTACGGTACGGCGTGGGCAGACAAAAAAGCACTTAACGCTTACCTGCAGCGCCTGGAAGAAGCCGCGAAACGCGACCACCGTAAAATCGGTAAACAGCTCGACCTGTACCATATGCAGGAAGAAGCGCCGGGTATGGTATTCTGGCACAACGACGGCTGGACCATCTTCCGTGAACTGGAAGTGTTTGTTCGTTCTAAACTGAAAGAGTACCAGTATCAGGAAGTTAAAGGTCCGTTCATGATGGACCGTGTCCTGTGGGAAAAAACCGGTCACTGGGACAACTACAAAGATGCAATGTTCACCACATCTTCTGAGAACCGTGAATACTGCATTAAGCCGATGAACTGCCCGGGTCACGTACAAATTTTCAACCAGGGGCTGAAGTCTTATCGCGATCTGCCGCTGCGTATGGCCGAGTTTGGTAGCTGCCACCGTAACGAGCCGTCAGGTTCGCTGCATGGCCTGATGCGCGTGCGTGGATTTACCCAGGATGACGCGCATATCTTCTGTACTGAAGAACAAATTCGCGATGAAGTTAACGGATGTATCCGTTTAGTCTATGATATGTACAGCACTTTTGGCTTCGAGAAGATCGTCGTCAAACTCTCCACTCGTCCTGAAAAACGTATTGGCAGCGACGAAATGTGGGATCGTGCTGAGGCGGACCTGGCGGTTGCGCTGGAAGAAAACAACATCCCGTTTGAATATCAACTGGGTGAAGGCGCTTTCTACGGTCCGAAAATTGAATTTACCCTGTATGACTGCCTCGATCGTGCATGGCAGTGCGGTACAGTACAGCTGGACTTCTCTTTGCCGTCTCGTCTGAGCGCTTCTTATGTAGGCGAAGACAATGAACGTAAAGTACCGGTAATGATTCACCGCGCAATTCTGGGGTCGATGGAACGTTTCATCGGTATCCTGACCGAAGAGTTCGCTGGTTTCTTCCCGACCTGGCTTGCGCCGGTTCAGGTTGTTATCATGAATATTACCGATTCACAGTCTGAATACGTTAACGAATTGACGCAAAAACTATCAAATGCGGGCATTCGTGTTAAAGCAGACTTGAGAAATGAGAAGATTGGCTTTAAAATCCGCGAGCACACTTTGCGTCGCGTCCCATATATGCTGGTCTGTGGTGATAAAGAGGTGGAATCAGGCAAAGTTGCCGTTCGCACCCGCCGTGGTAAAGACCTGGGAAGCATGGACGTAAATGAAGTGATCGAGAAGCTGCAACAAGAGATTCGCAGCCGCAGTCTTAAACAATTGGAGGAATAA

またAE000512側が

    gene            complement(763605..765527)
                    /locus_tag="TM_0740"
                    /old_locus_tag="TM0740"
    CDS             complement(763605..765527)
                    /locus_tag="TM_0740"
                    /old_locus_tag="TM0740"
                    /note="similar to GB:M36593 SP:P18256 PID:143764
                    PID:1565235 GB:AL009126 percent identity: 73.50;
                    identified by sequence similarity; putative"
                    /codon_start=1
                    /transl_table=11
                    /product="threonyl-tRNA synthetase"
                    /protein_id="AAD35821.1"
                    /translation="MKIKVKLPDGKEKEYDRGITPAEIAKELGIKKAIGAVVNGELWD
                    LKRPIENDCELRLVTLEDPEAPEFYRHTMAHILAQAVMRLYGKENVKLGIGPTIENGF
                    YYDFDIKNGRLTEEDLPKIEQEMKKIIKENLPIERKEISKEEARELFRDQPYKLELIE
                    EIEGNRVTIYRQGEFVDLCRGPHLPSTGIVKHFKLLSVSGAYWRGSEKNPMLTRVYGT
                    AFAKKEDLDNYLKFLEEAQRRDHRKLGPQLELFMLNTEYAPGMPFFLPKGVVVLNELM
                    KFSRELHRERGYQEIFTPLIMNEQLWKISGHWDHYAENMYFIEKDEERYAVKPMNCPG
                    HILVYKSRTVSYRDLPLRFFEFGRVHRYERSGVLHGLMRVRSFTQDDAHIFCTPDQIE
                    EEILGVLDLINTIYSQFGFTYRVELSTMPEDHMGDEAIWEKATTALKNALERAGLSYK
                    VNEGEGAFYGPKIDFHIRDSIGREWQCATIQLDFMMPEKFNVTYIGPDNKEHRAVMIH
                    RAIYGSLERFFGILIEHFAGAFPTWLAPIQVAVIPISEKHNDGAEKVARRISQEGFRV
                    FFDNRRETLGYRIRQAQTQKIPYMIIIGDKELESGKISVRTRTGKEIKDVDPEHFVET
                    LRNEVLSRKLELSMEG"

配列が

ATGAAGATAAAGGTGAAGCTTCCAGACGGAAAAGAAAAAGAGTACGACAGAGGCATCACACCGGCGGAGATCGCAAAAGAACTGGGCATAAAAAAGGCGATCGGAGCGGTTGTTAACGGTGAACTCTGGGACCTGAAAAGACCCATCGAAAATGACTGCGAACTGCGGCTTGTAACTCTGGAGGATCCCGAGGCACCCGAGTTCTACAGACACACCATGGCACACATCCTCGCACAGGCTGTGATGAGACTCTATGGAAAAGAAAACGTAAAACTGGGAATCGGCCCCACCATTGAGAACGGTTTTTACTACGATTTCGATATAAAGAACGGAAGGCTCACAGAGGAAGATCTCCCAAAGATAGAACAGGAGATGAAGAAGATAATAAAGGAAAATCTTCCTATAGAAAGAAAGGAGATCAGCAAAGAAGAGGCAAGAGAACTTTTCAGAGATCAGCCTTACAAACTCGAGCTAATAGAGGAGATCGAAGGCAACAGAGTAACGATCTACCGTCAGGGAGAGTTCGTGGATCTCTGCAGAGGTCCTCATCTTCCTTCGACAGGAATAGTGAAGCACTTCAAACTCCTCTCCGTTTCAGGCGCGTACTGGCGAGGCAGTGAGAAGAATCCCATGCTCACAAGAGTTTACGGAACGGCCTTCGCCAAAAAGGAAGACCTCGACAATTACTTGAAATTCCTTGAAGAAGCCCAGCGAAGGGACCACAGAAAATTGGGACCTCAACTCGAACTTTTCATGCTGAACACCGAGTACGCTCCCGGTATGCCGTTCTTCCTTCCAAAGGGAGTCGTAGTGCTCAACGAATTGATGAAATTTTCCAGAGAACTCCATCGTGAGAGAGGATACCAGGAGATCTTCACACCGCTCATAATGAACGAGCAGCTCTGGAAGATCTCCGGTCACTGGGATCACTACGCTGAGAACATGTACTTCATAGAGAAAGACGAGGAAAGGTACGCGGTGAAACCCATGAACTGTCCCGGTCACATTCTCGTGTACAAGAGCAGGACGGTTTCCTACAGGGACCTACCGCTGAGATTCTTCGAGTTTGGCCGGGTTCACAGATACGAGAGAAGCGGAGTTCTCCACGGTCTCATGAGGGTGAGATCCTTCACACAGGACGATGCACACATATTCTGCACTCCCGATCAGATCGAGGAAGAGATACTCGGCGTTCTCGATCTCATAAACACCATCTACAGCCAATTCGGCTTCACCTACAGAGTGGAACTCAGCACAATGCCCGAAGATCACATGGGTGACGAAGCGATTTGGGAGAAGGCAACAACCGCCTTGAAGAATGCCTTAGAAAGGGCCGGGCTCTCGTACAAAGTGAACGAGGGTGAAGGCGCTTTCTACGGTCCGAAAATAGACTTCCACATAAGAGATTCGATAGGAAGAGAGTGGCAGTGTGCCACGATTCAGCTGGATTTCATGATGCCCGAAAAGTTCAACGTCACTTACATAGGCCCAGACAACAAAGAGCACAGAGCCGTGATGATACACAGAGCGATATACGGTTCGCTCGAAAGGTTCTTTGGCATCCTGATAGAGCACTTCGCTGGAGCGTTTCCAACCTGGCTTGCACCGATCCAGGTGGCAGTTATTCCCATATCAGAAAAACACAACGACGGTGCGGAAAAGGTTGCAAGGAGGATCTCTCAGGAAGGGTTCAGGGTCTTCTTCGATAACCGCCGTGAAACTCTTGGATACCGTATCAGACAGGCTCAAACCCAGAAAATACCGTACATGATAATAATAGGTGACAAAGAGCTGGAGAGCGGGAAGATCTCCGTGAGAACAAGAACGGGAAAAGAGATAAAAGACGTAGATCCAGAACACTTCGTCGAAACTCTGAGAAACGAAGTTCTGAGCAGAAAGCTCGAACTTTCAATGGAGGGATAA

であった。

念のため、部分一致が正しく検出できているかの確認のため、次のようなテストを行った。

AE000512をデータベースにして、同じAE000512.fastaから一部を抜き出して更にシーケンスの前後を適当に切り捨てて短くした（下記TM_0005は長さ100を残し、0006は長さ200を残した）ものを用意した。これにより配列の途中部分が完全一致している例を作ったことになる。これをblastで比較した。

>ThermoTM_0005_Mod DNA helicase, putative
GAAGATAGACACAGAAATCGGTGTGGGTGATCTTGTTCTCATAAGCAAAGGAAACCCCCTCAAGAGTGATTACACGGGAACGGTTGTGGAGAAAGGAGAG
>ThermoTM_0006_Mod muconate cycloisomerase
GGAAGACATCGAAGCCGTGGAAGAGATAGCGAAAGTCACCCGTGGAGCGAAGTACATAGTTGACGCGAACATGGGATACACCCAGAAGGAGGCGGTGGAGTTCGCGAGAGCAGTATATCAAAAAGGAATAGACATCGCTGTGTACGAACAACCTGTGAGGAGGGAAGACATAGAAGGCTTGAAGTTCGTGAGGTTCCACT

このファイルAEshort.fastaに対して

blastn -db AE000512 -outfmt "10 qseqid sseqid pident length mismatch evalue bitscore" -evalue 100 < AEshort.fasta

を実行したところ、結果は

ThermoTM_0005_Mod,ThermoTM_0005,100.000,100,0,1.64e-48,185
ThermoTM_0006_Mod,ThermoTM_0006,100.000,200,0,9.45e-104,370

となり、正しく一致検出しただけでなく、一致部分の一致率100%で一致長さも100,200と正しく得られている。

となると、AP012030とAE000512の間での一致比較が１つしか出なかったことは、かなり確実に、部分一致するものが無かったと言えるのだろう。

↑

アミノ酸列に変換してから比較する †

GenBankアノテーションファイルから、CDSを抽出し、それぞれのCDSの＜アミノ酸の＞ FASTAファイルを作る。基本的には上と同じだが、配列を塩基配列からアミノ酸配列に 変換してから出力する。変換はBioPythonのTranslate機能。

#Read CDS and write Fasta -- AminoAcid version
import os, sys
import numpy as np
import pandas as pd
from Bio import SeqIO
from Bio.Seq import *
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import SeqFeature, FeatureLocation
from collections import OrderedDict
import pickle
from ReadCDSforComparison import ReadCDS
def OutputCDS(tab1, tab2):
    with open('test-Clustalo.fasta', 'w') as cf:
        for i,r in tab1.iterrows():
            print('>Ecoli_%s %s\n%s\n' % (r['locus_tag'].strip(), r['product'], r['seq']), file=cf)
        for i,r in tab2.iterrows():
            print('>Thermo_%s %s\n%s\n' % (r['locus_tag'].strip(), r['product'], r['seq']), file=cf)
    return

def OutputCDSforBlastp(tab1, tab2):
    with open('AP012030CDSAA.fasta', 'w') as cf:
        for i,r in tab1.iterrows():
            AAseq = r['seq'].translate(table=11, stop_symbol="")
            print('>%s %i %s\n%s\n' % (r['locus_tag'].strip(), len(AAseq), r['product'], AAseq), file=cf)

    with open('AE000512CDSAA.fasta', 'w') as cf:
        for i,r in tab2.iterrows():
            AAseq = r['seq'].translate(table=11, stop_symbol="")
            print('>Th_%s %i %s\n%s\n' % (r['locus_tag'].strip(), len(AAseq), r['product'], AAseq), file=cf)

    return
    
if __name__ == '__main__':
    ReviseEcoliPickle = False
    #ReviseEcoliPickle = True
    ReviseThermoPickle = False
    #ReviseThermoPickle = True
    
    Ecoligbfile = 'AP012030.gb'
    EcoliCDSfile = 'AP012030CDS.pickle'
    if os.path.exists(EcoliCDSfile) and not ReviseEcoliPickle:
        with open(EcoliCDSfile, 'rb') as pf:
            EcoliCDS = pickle.load(pf)
    else:
        EcoliCDS = ReadCDS(Ecoligbfile)
        with open(EcoliCDSfile, 'wb') as pwf:
            pickle.dump(EcoliCDS, pwf)

    Thermogbfile = 'AE000512.gb'
    ThermoCDSfile = 'AE000512CDS.pickle'
    if os.path.exists(ThermoCDSfile) and not ReviseThermoPickle:
        with open(ThermoCDSfile, 'rb') as pf:
            ThermoCDS = pickle.load(pf)
    else:
        ThermoCDS = ReadCDS(Thermogbfile)
        with open(ThermoCDSfile, 'wb') as pwf:
            pickle.dump(ThermoCDS, pwf)

    OutputCDSforBlastp(EcoliCDS, ThermoCDS)
   
    print('complete')

これで得られるAP012030CDSAA.fastaとAE000512CDSAA.fastaを使ってblastpを動作させる。 この手順は上記のDNA配列時と同じ。

まず、配列データベースを作る。

makeblastdb -in AE000512CDSAA.fasta -parse_seqids -blastdb_version 5 -title "ThermoAA" -out AE000512AA -dbtype prot

結果は

Building a new DB, current time: 09/05/2019 06:51:46
New DB name:   /home/yamanouc/src/Kishimoto2017/test-Blast/AE000512AA
New DB title:  ThermoAA
Sequence type: Protein
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 1846 sequences in 0.10688 seconds.

このデータベースを使って、blastp を起動し、AP012030側のCDSAAのfastaファイルAP012030CDSAA.fastaを入力として与える。今回はかなり多数のヒットが見られるので、E-valueの閾値はデフォルト通りE<10とする。

blastp -db AE000512AA -outfmt "10 qseqid sseqid pident length mismatch evalue bitscore" -out AP012030AAblast.out -num_threads 32 -evalue 10 < AP012030CDSAA.fasta

逆方向で試す。データベースを作る。

makeblastdb -in AP012030CDSAA.fasta -parse_seqids -blastdb_version 5 -title "EcoliAA" -out AP012030AA -dbtype prot

AE000512CDSAA.fastaと比較する。

blastp -db AP012030AA -outfmt "10 qseqid sseqid pident length mismatch evalue bitscore" -out AE000512AAblast.out -num_threads 32 -evalue 10 < AE000512CDSAA.fasta

結果は、多数の類似シーケンスが抽出された。 |データベースvs対象 | E<10 のマッチング件数|

AE000512 vs AP01203047949
AP012030 vs AE000512147190

かなり多いので、さらに絞り込む必要がある。

blastpの起動は、biopythonを介して行うことができる。これのメリットは、入出力の接続が（ファイル経由に比べて）やや便利かな、というところ。さほど便利でもないか？

とりあえず、次のようにするとできる。

from Bio.Blast.Applications import NcbiblastpCommandline
#help(NcbiblastpCommandline)

blastp_cline = NcbiblastpCommandline(db='AP012030AA', query='AE000512CDSAA.fasta', 
                                     outfmt='"10 qseqid sseqid pident length mismatch evalue bitscore"',
                                     out='AE000512AAblast2.out2', num_threads=32, evalue=0.001)
blastp_cline   # これで実行？

print(blastp_cline)
# これによって以下のように書きだす。
blastp -out AE000512AAblast2.out2 -outfmt "10 qseqid sseqid pident length mismatch evalue bitscore" -query AE000512CDSAA.fasta -db AP012030AA -evalue 0.001 -num_threads 32
#stdout, stderr = blastp_cline()   # これは要らない？
#print(stdout)
print('complete')

これによって、-outで指定したファイル AE000512AAblast2.out2 が生成される。

なお、コマンドラインでは bitscore による閾値絞り込みはできない。evalueのみ。したがって、たとえば -evalue=0.001 とかしておいて、出力ファイルの上で絞り込む必要がある。

↑

結果の絞り込み †

blastpで得られたそれぞれの一致候補に対して、GenBankアノテーションから LocusIDを頼りにして、product情報と配列の長さ（アミノ酸列の長さ）情報を追記する。

#Read blast-out and augment with length(AA) and product
import os, sys
import numpy as np
import pandas as pd
from collections import OrderedDict
import pickle

def read_blast_out(blast_out_file):
    if not os.path.exists(blast_out_file):
        print('blast-out file', blast_out_file, 'does not exist')
        return
    df = pd.read_csv(blast_out_file, sep=',', header=None)
    df.columns = ['locus_tag_xx', 'locus_tag_yy', 'pident', 'length', 'mismatch', 'evalue', 'bitscore']
    return(df)

if __name__ == '__main__':
    EcoliCDSfile = 'AP012030CDS.pickle'
    if os.path.exists(EcoliCDSfile):
        with open(EcoliCDSfile, 'rb') as pf:
            EcoliCDS = pickle.load(pf)
    else:
        print('EcoliCDS file not found', EcoliCDSfile)
        return
    
    ThermoCDSfile = 'AE000512CDS.pickle'
    if os.path.exists(ThermoCDSfile):
        with open(ThermoCDSfile, 'rb') as pf:
            ThermoCDS = pickle.load(pf)
    else:
        print('ThermoCDS file not found', ThermoCDSfile)
        return

    # Read blast-out file
    APDF = read_blast_out('AP012030AAblast.out')
    print('----')
    print(EcoliCDS.head())
    print('----')
    print(APDF.head())
    print('----')
    ThermoCDS['locus_tag'] = ThermoCDS['locus_tag'].map(lambda u: 'Th_'+u)
    ThermoCDS = ThermoCDS.drop(columns='AAseq')
    print(ThermoCDS.head())
    print('----')
    AugAPDF = pd.merge(APDF, EcoliCDS, left_on='locus_tag_xx', right_on='locus_tag')
    AugAPDF = pd.merge(AugAPDF, ThermoCDS, left_on='locus_tag_yy', right_on='locus_tag')
    print(AugAPDF.head())
    AugAPDF = AugAPDF.sort_values(['locus_tag_yy', 'evalue', 'locus_tag_xx'])
    AugAPDF.to_excel('AP012030AAblast.xlsx')
    ###
    AEDF = read_blast_out('AP012030AAblast.out')
    AugAEDF = pd.merge(AEDF, ThermoCDS, left_on='locus_tag_yy', right_on='locus_tag')
    AugAEDF = pd.merge(AugAEDF, EcoliCDS, left_on='locus_tag_xx', right_on='locus_tag')
    print(AugAEDF.head())
    AugAEDF = AugAEDF.sort_values(['locus_tag_xx', 'evalue', 'locus_tag_yy'])
    AugAEDF.to_excel('AE000512AAblast.xlsx') 

    print('complete')

結果の先頭部分はこんな（添付ファイル）感じになった

↑

全体をまとめた処理を考える †

１）AP012030AAのデータベースは予め作っておく

２）対象好熱菌gbから、CDS部分のアミノ酸列を取出して、FASTAファイルに作る。
afile.fasta <- getAA(afile.gb)

３）blastpで、db=AP012030AAとして query=afile.fasta として、出力を out=afile.blt

４）blastpの出力から、bitscore>=200 のエントリーだけをフィルタしてafile200.bltを出力する。

５）それに、両方のデータのproductを付加する

６）さまざまな対象好熱菌の出力を集める。

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2019-10-22 (火) 17:05:04 (1254d)