search

soedinglab / MMseqs2

MMseqs2: ultra fast and sensitive search and clustering suite
https://mmseqs.com
GNU General Public License v3.0
1.31k stars 184 forks source link

Incorrect identity reported in m8 file #116

Open josemduarte opened 5 years ago

josemduarte commented 5 years ago

Expected Behavior

A sequence search that match itself in target db should report 100% identity

Current Behavior

Identity reported is 91.4%:

P30340_1-122    UniRef90_P30340 0.914   122     10      0       1       122     1       122     3.11E-76        235

Steps to Reproduce (for bugs)

My test_target_db.fa is:

>UniRef90_A9WC40 Succinyl-CoA--L-malate CoA-transferase alpha subunit n=3 Tax=Chloroflexus TaxID=1107 RepID=SMTA_CHLAA
MAKASRLTRSTGQPTEVSEGQVTGTSEMPPTGEEPSGHAESKPPASDPMSTPGTGQEQLP
LSGIRVIDVGNFLAGPYAASILGEFGAEVLKIEHPLGGDPMRRFGTATARHDATLAWLSE
ARNRKSVTIDLRQQEGVALFLKLVAKSDILIENFRPGTMEEWGLSWPVLQATNPGLIMLR
VSGYGQTGPYRRRSGFAHIAHAFSGLSYLAGFPGETPVLPGTAPLGDYIASLFGAIGILI
ALRHKEQTGRGQLIDVGIYEAVFRILDEIAPAYGLFGKIREREGAGSFIAVPHGHFRSKD
GKWVAIACTTDKMFERLAEAMERPELASPELYGDQRKRLAARDIVNQITIEWVGSLTRDE
VMRRCLEKEVPVGPLNSIADMFNDEHFLARGNFACIEAEGIGEVVVPNVIPRLSETPGRV
TNLGPPLGNATYEVLRELLDISAEEIKRLRSRKII
>UniRef90_A9WC39 Succinyl-CoA--L-malate CoA-transferase beta subunit n=7 Tax=Chloroflexus TaxID=1107 RepID=SMTB_CHLAA
MDGTTTTLPLAGIRVIDAATVIAAPFCATLLGEFGADVLKVEHPIGGDALRRFGTPTARG
DTLTWLSESRNKRSVTLNLQHPEGARVFKELIAHSDVLCENFRPGTLEKWGLGWDVLSKI
NPRLIMLRVTGYGQTGPYRDRPGFARIAHAVGGIAYLAGMPKGTPVTPGSTTLADYMTGL
YGCIGVLLALRHREQTGRGQYIDAALYESVFRCSDELVPAYGMYRKVRERHGSHYNEFAC
PHGHFQTKDGKWVAISCATDKLFARLANAMGRPELASSSVYGDQKVRLAHASDVNEIVRD
WCSSLTRAEVLERCYATATPAAPLNDIADFFGDRHVHARRNLVAIDAEDLGETLIMPNVV
PKLSETPGSIRSLGPKLGEHTEEVLKEILGMCDEQINDLRSKRVI
>UniRef90_P9WMI4 HTH-type transcriptional repressor SmtB n=18 Tax=Mycobacterium tuberculosis complex TaxID=77643 RepID=SMTB_MYCTO
MVTSPSTPTAAHEDVGADEVGGHQHPADRFAECPTFPAPPPREILDAAGELLRALAAPVR
IAIVLQLRESQRCVHELVDALHVPQPLVSQHLKILKAAGVVTGERSGREVLYRLADHHLA
HIVLDAVAHAGEDAI
>UniRef90_P30340 Transcriptional repressor SmtB n=5 Tax=Bacteria TaxID=2 RepID=SMTB_SYNE7
MTKPVLQDGETVVCQGTHAAIASELQAIAPEVAQSLAEFFAVLADPNRLRLLSLLARSEL
CVGDLAQAIGVSESAVSHQLRSLRNLRLVSYRKQGRHVYYQLQDHHIVALYQNALDHLQE
CR
>UniRef90_H2E7T8 Sterol methyltransferase-like 1 n=1 Tax=Botryococcus braunii TaxID=38881 RepID=SMTL1_BOTBR
MASELFATYYPRVVEAAQQLAPWQIAAGVTAAVVIGGYIWIITELRSPRRTGTSLFKLSG
GGIKKHDVAKFMDGYEKSYKTQEDGALTWHHISKEDSVKMVNTFYDLVTDAYEWAWDISF
HFSCRPVWANFAQAQVLHECRIANLARIQPGMKVIDVGTGVGNPGRTIASLTGAHVTGVT
INAYQIKRALHHTKKAGLLDMYKPVQADFTDMPFADESFDAAFAIEATCHAPKLEQVYAE
VYRVLKPGAYFAVYEAVSKPNFDPKNKRHVEIINSLVYGNGIPDMRTWKEAEEAGKKVGF
KLHFSYDAGEASSVLAPWWERPRNLVNTGVIAYTKFAIKVCDKIGILPRDYAKFAKCVGD
CIPDAVESGELGIFTPMYVYVWQKPEKST
>UniRef90_A8MU46 Smoothelin-like protein 1 n=7 Tax=Simiiformes TaxID=314293 RepID=SMTL1_HUMAN
MEQKEGKLSEDGTTVSPAADNPEMSGGGAPAEETKGTAGKAINEGPPTESGKQEKAPAED
GMSAELQGEANGLDEVKVESQREAGGKEDAEAELKKEDGEKEETTVGSQEMTGRKEETKS
EPKEAEEKESTLASEKQKAEEKEAKPESGQKADANDRDKPEPKATVEEEDAKTASQEETV
VEDEAKAEPKEPDGKEEAKHGAKEEADEPGSPSEEQEQDVEKEPEGGAGVIPSSPEEWPE
SPTGEGHNLSTDGLGPDCVASGQTSPSASESSPSDVPQSPPESPSSGEKKEKAPERRVSA
PARPRGPRAQNRKAIVDKFGGAASGPTALFRNTKAAGAAIGGVKNMLLEWCRAMTKKYEH
VDIQNFSSSWSSGMAFCALIHKFFPDAFDYAELDPAKRRHNFTLAFSTAEKLADCAQLLD
VDDMVRLAVPDSKCVYTYIQELYRSLVQKGLVKTKKK

and my test_query.fa is one of the sequences in the target:

>P30340_1-122
MTKPVLQDGETVVCQGTHAAIASELQAIAPEVAQSLAEFFAVLADPNRLRLLSLLARSEL
CVGDLAQAIGVSESAVSHQLRSLRNLRLVSYRKQGRHVYYQLQDHHIVALYQNALDHLQE
CR

The commands I ran are:

/opt/mmseqs2/bin/mmseqs createdb test_target_db.fa test_target_dbDB
/opt/mmseqs2/bin/mmseqs createdb test_query.fa test_queryDB

/opt/mmseqs2/bin/mmseqs search test_queryDB test_target_dbDB test_resultDB tmp/
/opt/mmseqs2/bin/mmseqs convertalis test_queryDB test_target_dbDB test_resultDB test_result.m8

MMseqs Output (for bugs)

The .m8 output is:

P30340_1-122    UniRef90_P30340 0.914   122     10      0       1       122     1       122     3.11E-76        235
P30340_1-122    UniRef90_P9WMI4 0.326   90      60      0       31      120     43      132     1.71E-11        49

Your Environment

martin-steinegger commented 5 years ago

Thank you for reporting this. We estimate the sequence identity in default by the score per column. This is useful for lower identities. If you want the real sequence identity please use -a. See also https://github.com/soedinglab/MMseqs2/wiki#how-does-mmseqs2-compute-the-sequence-identity

josemduarte commented 5 years ago

Thanks @martin-steinegger ! Sorry I missed that in the docs. A possible usability improvement could be not outputting identity values in -a false case.