Thanks for help, but still, I don't get the same result as blast. Okay, I'll give you the example on which I'm trying to do it in the hope you can get something out of it. I run blastp on Uniprot, with this query: Query= d1r5la1 a.5.3.1 (A:25-90) Alpha-tocopherol transfer protein {Human (Homo sapiens) [TaxId: 9606]} Length=66 I took that query from SCOP/Astral database. The first significant alignment I get from blast is for > sp|P49638|TTPA_HUMAN. and blast says this: Score = 136 bits (343) and Expect = 1e-39. Okay. When I plug in the stuff you told me to, I get: 1.448e-31. That is nowhere close the blast's value. Assuming database length is 191,240,774, K is 0.041, and lambda is 0.267. (Although I didn't need them right now.)

I think that recent versions of blastp use a slightly different stats model than the classical Altschul one. I faced similar discrenpancies between e-values from blastp and those I computed myself the same way you did.

What is your blast version ? If it is quite recent, you should try version 2.2.26 or 2.2.27 for instance (I guess they still use the old model). Otherwise, you could set the environment variable OLD_FSC to 1 before launching your blast command; it will switch back to the old model.