BLAST revisit

上星期(2014/2/26)又到醫學院講述序列比對的重要觀念與操作細節
回想起去年上完課後寫下的文章，把它翻出來與大家分享。

今年，想藉機會進一步討論BLAST的e-values。

讓我們先用以下序列(query)做個試驗：

>FASTA format

MAHTTSLLSGASSAAASPLSSRQYSVASSPRLSGDIGRGLFAIVVLLLRMSSVYGVIDRL VVQTSSGPVRGRSVTVQGREVHVYTGIPYAKPPLDDLRFRKPVPAEPWHGVLDATRLPAT CVQERYEYFPGFSGEEIWNPNTNVSEDCLYINVWAPAKARLRHGRGANGGEHSNKADTDH LIHNGNPQNTTNGLPVLIWIYGGGFMTGTATLDIYNADIMSAVGNVIVASFQYRVGAFGF LHLSPAMPGYEEEAPGNVGLWDQALAIRWLKTNAHAFGGNPEWMTLFGESAGSSSVNAQL VSPVTAGLVKRGMMQSGTMNAPWSHMTSEKAVEIGKALINDCNCNASLLSENPQAVMACM RAVDAKTISVQQWNSYSGILSFPSAPTIDGAFLPDHPMKMMETADLRGYDILMGNVRDEG TYFLLYDFIDYFDKDEATSLPRDKYLEIMNNIFGKVTQAEREAIIFRHTSWVGNPGLENQ QQIGRAVGDHFFTCPTNEYAQALAERGASVHYYYFTHRTSTSLWGEWMGVLHGDEIEYFF GQPLNTSLQYRQVERELGKRMLNAVIEFAKTGNPATDGEEWPNFTKKDPVYYVFSTDDKE EKLQRGPLEGRCAFWNEYLREVRKWGSQCELKPSSASSLQQKQQHLLLQQRSIVTFMLTL

SLVLGIPSVNAFF 

1. 首先選擇 blastp + nr (13,279,552,322 bps)
2. 再選擇 blastp + swiss-prot (170,826,989 bps)
3. 再選擇 blastp + pdb (16,752,417 bps)

以搜尋結果中的human Acylcholine acylhydrolase (Swiss-prot ID: P06276, PDB ID: 2WIF:A)為例：

Score	E-value	Ident	Accession	Database
397	8e-127	39%	2WIF_A (length: 529)	nr
397	2e-127	38%	P06276.1 (length: 602)	swiss-prot
397	1e-129	39%	2WIF_A (length: 529)	pdb

可以看出，雖然alignment score都一樣(397)，但e-value卻有數十倍到數百倍的差異。其差異主要來自於database的大小(nr > swiss-prot > pdb)，一般而言，較小的database (e.g. pdb)會得到較小的e-value!

當然，如果看到上述這麼小的e-values，自然沒什麼好擔心的。但如果得到像0.01這樣的e-value，就要稍微注意一下。

接下再做另一個測試，假設我們只有局部的序列：

>FASTA format

CVQERYEYFPGFSGEEIWNPNTNVSEDCLYINVWAPAKARLRHGRGANGGEHSNKADTDH LIHNGNPQNTTNGLPVLIWIYGGGFMTGTATLDIYNADIMSAVGNVIVASFQYRVGAFGF

同樣用blastp對swiss-prot進行序列比對，這次對P06276的結果變成：

Score	E-value	Ident	Accession	Database
103	5e-25	45%	P06276.1	swiss-prot

這次，e-value變大的主因來自於score變小!