Glitches in GFF3 match format output (-O2)

[ohdongha]

Dear Osamu, @ogotoh

I followed your instruction from issue 77 to successfully map and align protein sequences to a genome. However, the "match format" GFF3 appears to have some glitches.

1. I downloaded the genome sequence from NCBI:

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/034/640/455/GCF_034640455.1_iyDiaLong2/GCF_034640455.1_iyDiaLong2_genomic.fna.gz

2. The protein sequence input (toy_example.fa) was as follows:

>gi|1983991430|ref|XP_039310866.1| ATP-binding cassette sub-family G member 1 [Solenopsis invicta]
MFRKIPKGSNKNNAARMPYDECRFDLQDAEKLLKSEKWNSRTYIEFNDLCYSIRNHKGSE
NTILHNVTGHFEPGKITVIIGPSGAGKTTLMKIISGKRSMDIKGTLTVNNDEWNKGMFRK
HVCYVPQQFDLLPYLTTRETLCIAARLKLDVNQNKQEINTVVNDIAENLKLSNCLDTLAN
RLSGGEKKRLSIGVQMITKSSVFLLDEPTSGLDSAASYQLLSILHNMAKANCTIVCAIHQ
PSSRMISLFDDIMVLNRGRCLYCGPKSEILSTFNIAGFTCPSFYNITEFVLDVITEQEDE
DLKNLQKISRDAYEKFRPRSEHKNKLNSLTDLKQNFERDNTITTINSKMQEKSTWQQQKI
LFSRTLIGIKRDSTILKIASHILSGLLFGSMFYNFGDDATKVQSNIFCLFITPVYLFFMH
SLQTVMTISSETAVFLQEYLNNWYSFRAYYSIKVLSDILMQVLCTTSYLLILYYMTGQPM
VFNRIIQTWTICLLIMILGQTYGIFVGTAFGAKIGVVLIVAINVPLMMFAGFYVKIAELP
KYIQPLGYISYFRYSFEGMM
>gi|1321349347|ref|XP_023316908.1| nitric oxide synthase-like [Trichogramma pretiosum]
MSTPPHGTEVRKSNEILEHAKDFLNQYFTSIKRVNSESHKSRWKNIEKEVLSTGSYQLTE
TELIFGAKLAWRNASRCIGRIQWSKLQVFDCRYVTTTSGMFEALCNHIKYSTNKGNIRSA
ITVFPQRTDGKHDFRVWNPQLISYAGYKNDDGSFIGDAGNVEFTEVCVKLGWKSARTQFD
VLPLVLSANGHDPDYFDIPEELILEIPLVHPSFEWFGSLGLKWYAVPAVSNMVFDCGGLE
FTAAPFNGWYMSTEIGSRDLCDVNRYNLLEVVAQKMDLDTKVTSSLWKDKAMIEVNLAVL
YSFQIKNVTIVDHHTASDSFMKHYENEIRLRHGCPADWVWIIPPLSGSATAVFHQEMALY
HLKPSYNTQEPAWKTHIWRKERNKYSTTKQLRRKFHFKQIARAVKFTSKLFGRALSRRIK
ATVLYATETGKSQAYADKLGELFGHAFHSQVLCMADYDITKIEHEALLLVITSTFGNGDP
PENGEVFAQNLYEMKIDDPFITSDTQFDIRTSKSFIKVNSQSGITKSKKINKLTSLKELD
YKADENFGPLSNVR
>gi|2322517383|ref|XP_051166906.1| lebercilin-like protein isoform X1 [Leptopilina boulardi]
MQTSKTPLPSLCLSSRLQREENYVPGIKSDKKEHSKSAAFLLSNRKQLNSLNGGPYNSRN
RCYVQPIHLKQSNKSAFIQKVVSAKMLQIKALQNELTDAHYHLNELANENRLLKTLQTRQ
DSALRKYEGTNAELPRIINTHHEELRVLQTKYKKLKSQYRSNCELLKEKDVKIHIIQDQN
KHLLKLSKDRNLEERESLQRKTSDLIFKIEQQDKTINTLNRKLAIETKHLKHQLHLEVIK
HKETQKKLQETFNKLKSTEELLENKEKKLLHGKHFNLIRTMQNIKSTSHPNLDNNRLGNS
DGDNFNTPKVSLNNLENKKDSLIIEDVEKFIQQLNLSKISNRMDSIATISQCHITEIIEH
NKDNNYLTQESFERSNENKIVISKQDNGMDNNSYNMNDNYVNNEIKDSEISLITTTESVS
SLQNNDSNSPRKEKLHLKYSISEESEPEDFADELNDEIKSEICNSVKTTNNTTETLVFQK
KNNYDKERLLAAMKAIDDNENIE
>gi|2322533789|ref|XP_051175486.1| GPI transamidase component PIG-S [Leptopilina boulardi]
MSNSQEDAVGHDTATDEKYRTYASISFAVLLLGVGIPLWWYTTAVPRVPLPYSGIAALSH
HETKIQTDIIVAALSKERAEQLAQEIKSFFEDSEIYSLRVRFEVVSRNLMSSAFTTNELE
TIAYSFNLKMGNLLLLEALNLNEEILVGSNRTVFFAPDTDVKKVTEIVSKWILRDKSLAL
TKNALTHPTEYSLDKENRRRFPTSPAYDVLLTLVNPDPEKLKVNWSLPTIADVYLQPFLN
EISLVADFHIKSQWLYLVPLNVKPKEVYDSSILGRHFALSEDILPQLITPLEKKLASQVS
LHPCINLVAYVVPCEEGPLHIYSRNGHRAQTVNIEAFHSPRWGGVVIINPPVDACISTKP
VEITPNDATVIGVFFAQLRLLLGIPDPTPIPGAAVVPLFGFKPREWELDTLLRVRAVEQL
TSAKLTLQSLAQLLEEISNIVITETVGNRIKNALNFVEESSALLQSGNLRDGFLKSKEAF
ATAEAAFTDPSLLALLYFPDDQKYAVYIPLFLPVMIPVLLSLKNIRQYYFPR
>gi|645029614|ref|XP_008217361.1| transforming acidic coiled-coil-containing protein 3 isoform X1 [Nasonia vitripennis]
MDFLNRLLNRPTPKPAICHGSPTRQQHESAGEAALPSSRPQSDSDHNLIRFSRDFSEAAS
NSSGPSSLVPESLSSASSFQSLASTLSSKSSSGICADSAPGSPDKVNSGNGFDHSYVNVD
IPDLDDLISACANINLNNATVLGGLLPGENRAGDESPDSSGGDLSQDQTFASAADDTLAD
ETQVLVEETIQTFKTCPPLEVTQVTSAETHVFKAPISVDEQSSVIPYLSFEESVLRKHQE
FVDDGAEFKSALTSFQATQSITVSSTVDTRRGSYSPTPCEEIDVAQVAKITCKSPKDLDI
AQESLIVVVDEVVEKKEQIVEQDFILGTPDASFEEVQAFNSTIILNTDLASEEPTVLSNI
KNEVVDKEEKKPVISLDEVPVGVKTKGSEIVQETIPNILNHTHDNAQDTGILNKTQDIVQ
DANVLNNTQDIVQNTNILNQTQDIVQEKTILNKTQDIVQEESNLIKTELIAQEVNKSVNK
NEVIIRQISDLSKTESVQSKTQHFTKEDHFVEKKEILKHNIGSYNEIDERPIKPQIQEAQ
SPKIIESNEENTCRTPQKLEETFPKPEAKSPKQLETFANEEAQNLNLINTEAVASTPLPV
KNLNKAKTPKTPKTKSALFTNEPKSPKPPAQTSKTPIKSEEIKEEKNKLNFDQTLAINKS
LEQSSFSATELITNLETTSDVLHENLCIAPIQEPEYEAFKPIQQSTTLPFPDFKKLQEAA
LSVAKEIDESLEKIEETEQFVSATTELFGDPSALDFLTQVGSQHPTKRGERYDSLYMKFD
PLAQRLSMLPQKNLPTAPSVDEEEEKEKEAKVNKPSSEMGTPKKNPALAAIDRLLFYSPC
PPSATESAVAEAKQKKEEAEKQESPVVPFVDEKMAKELELVRTTVLQLEAELEKQKQEHR
EQLEKQKAAYEEKVTHMQTKFTQMQSQLTQEVKCKGQMTVVVEEYEKSISRLIAEREKDR
TSFEVEKAKLQEELQVANQHLSNTEAAFSDVHSKYERLKNVVSAYKSNETVLKESIAENQ
ETIKTLENRYEQLKSHAMAQLEKANLELSAIRKQNEAETVKLKALVRKAELKSKSLAETV
EQKIKENKELTQIIDEMIARVD
>gi|2322488769|ref|XP_051176403.1| acireductone dioxygenase isoform X1 [Leptopilina boulardi]
MVNIWYMDNSNADQRFDHQKIPQHFLSTQELLELTGVEYFKINASNHKTDTVLRDLKKTR
NYTYEDEITCSKEYLENYEEKLKQFFVEHLHTDEEIRLVLSGSGYFDVRDKDDQWIRIEV
KEGDLIIIPSGIYHRFTLDMNNFIKAKRYFVGEPIWLPYKRPADEMDCRKQYLNRLKTGF

3. Codes to run Spaln (version 3.0.6a <240916>)

spaln -W gnm_test -KP GCF_034640455.1_iyDiaLong2_genomic.fna.gz
spaln -Q7 -d gnm_test -A0 -O2 -o toy_example.gff3 -t8 -po toy_example.fa 2> toy_example_Q7.stderr.log

Then the output toy_example.gff3 includes multiple places where the start and end locations appear swapped:

##gff-version   3
##sequence-region       NC_087242.1 2068481 2140160
NC_087242.1     ALN     nucleotide_to_protein_match     2103841 2103963 222     +       .       ID=match00001;Name=NC_087242.1_2106;Target=XP_051176403.1  1 41 +;Gap=M41 
NC_087242.1     ALN     nucleotide_to_protein_match     2104319 2104438 171     +       .       ID=match00001;Name=NC_087242.1_2106;Target=XP_051176403.1  42 81 +;Gap=M40 
NC_087242.1     ALN     nucleotide_to_protein_match     2105780 2105959 320     +       .       ID=match00001;Name=NC_087242.1_2106;Target=XP_051176403.1  82 141 +;Gap=M60 
NC_087242.1     ALN     nucleotide_to_protein_match     2108306 2108422 222     +       .       ID=match00001;Name=NC_087242.1_2106;Target=XP_051176403.1  142 180 +;Gap=M39 
NC_087242.1     ALN     nucleotide_to_protein_match     2108444 2108443 5       +       .       ID=match00001;Name=NC_087242.1_2106;Target=XP_051176403.1  181 180 +;Gap=
##sequence-region       NC_087229.1 2949121 3031040
NC_087229.1     ALN     nucleotide_to_protein_match     2982562 2983262 690     +       .       ID=match00002;Name=NC_087229.1_2983;Target=XP_008217361.1  803 1043 +;Gap=M17 I1 M28 I1 M16 I1 M33 D3 M15 I7 M122 
NC_087229.1     ALN     nucleotide_to_protein_match     2984212 2984390 244     +       .       ID=match00002;Name=NC_087229.1_2983;Target=XP_008217361.1  1044 1102 +;Gap=M59 
NC_087229.1     ALN     nucleotide_to_protein_match     2984409 2984407 -3      +       .       ID=match00002;Name=NC_087229.1_2983;Target=XP_008217361.1  1103 1102 +;Gap=
##sequence-region       NC_087227.1 8243201 8314880
NC_087227.1     ALN     nucleotide_to_protein_match     8278318 8278415 221     +       .       ID=match00003;Name=NC_087227.1_8280;Target=XP_023316908.1  1 33 +;Gap=M33 
NC_087227.1     ALN     nucleotide_to_protein_match     8278534 8278788 438     +       .       ID=match00003;Name=NC_087227.1_8280;Target=XP_023316908.1  34 118 +;Gap=M85 
NC_087227.1     ALN     nucleotide_to_protein_match     8278895 8279036 235     +       .       ID=match00003;Name=NC_087227.1_8280;Target=XP_023316908.1  119 165 +;Gap=M47 
NC_087227.1     ALN     nucleotide_to_protein_match     8279114 8279265 232     +       .       ID=match00003;Name=NC_087227.1_8280;Target=XP_023316908.1  166 212 +;Gap=M10 D4 M37 
NC_087227.1     ALN     nucleotide_to_protein_match     8279406 8279580 331     +       .       ID=match00003;Name=NC_087227.1_8280;Target=XP_023316908.1  213 270 +;Gap=M58 
NC_087227.1     ALN     nucleotide_to_protein_match     8279708 8279809 151     +       .       ID=match00003;Name=NC_087227.1_8280;Target=XP_023316908.1  271 304 +;Gap=M34 
NC_087227.1     ALN     nucleotide_to_protein_match     8279939 8280231 520     +       .       ID=match00003;Name=NC_087227.1_8280;Target=XP_023316908.1  305 402 +;Gap=M98 
NC_087227.1     ALN     nucleotide_to_protein_match     8280308 8280625 496     +       .       ID=match00003;Name=NC_087227.1_8280;Target=XP_023316908.1  403 508 +;Gap=M106 
NC_087227.1     ALN     nucleotide_to_protein_match     8280780 8280923 140     +       .       ID=match00003;Name=NC_087227.1_8280;Target=XP_023316908.1  509 554 +;Gap=M30 D2 M16 
NC_087227.1     ALN     nucleotide_to_protein_match     8282427 8282430 -30     +       .       ID=match00003;Name=NC_087227.1_8280;Target=XP_023316908.1  555 554 +;Gap=D1 
##sequence-region       NC_087230.1 3747841 3819520
NC_087230.1     ALN     nucleotide_to_protein_match     3797123 3797168 110     -       .       ID=match00004;Name=NC_087230.1_3788;Target=XP_051175486.1  1 16 +;Gap=M6 I1 M9 
NC_087230.1     ALN     nucleotide_to_protein_match     3783586 3783816 298     -       .       ID=match00004;Name=NC_087230.1_3788;Target=XP_051175486.1  17 93 +;Gap=M77 
NC_087230.1     ALN     nucleotide_to_protein_match     3783167 3783361 199     -       .       ID=match00004;Name=NC_087230.1_3788;Target=XP_051175486.1  94 159 +;Gap=M48 I1 M17 
NC_087230.1     ALN     nucleotide_to_protein_match     3781962 3782180 282     -       .       ID=match00004;Name=NC_087230.1_3788;Target=XP_051175486.1  160 232 +;Gap=M73 
NC_087230.1     ALN     nucleotide_to_protein_match     3781711 3781899 277     -       .       ID=match00004;Name=NC_087230.1_3788;Target=XP_051175486.1  233 295 +;Gap=M63 
NC_087230.1     ALN     nucleotide_to_protein_match     3780976 3781262 314     -       .       ID=match00004;Name=NC_087230.1_3788;Target=XP_051175486.1  296 387 +;Gap=M37 D1 M35 D3 M20 
NC_087230.1     ALN     nucleotide_to_protein_match     3780557 3780906 442     -       .       ID=match00004;Name=NC_087230.1_3788;Target=XP_051175486.1  388 503 +;Gap=M7 D1 M109 
NC_087230.1     ALN     nucleotide_to_protein_match     3779079 3779165 150     -       .       ID=match00004;Name=NC_087230.1_3788;Target=XP_051175486.1  504 532 +;Gap=M29 
NC_087230.1     ALN     nucleotide_to_protein_match     3779030 3779033 -38     -       .       ID=match00004;Name=NC_087230.1_3788;Target=XP_051175486.1  533 532 +;Gap=D1 
##sequence-region       NC_087233.1 6512641 6676480
NC_087233.1     ALN     nucleotide_to_protein_match     6528433 6528467 91      -       .       ID=match00005;Name=NC_087233.1_6523;Target=XP_051166906.1  1 15 +;Gap=M9 I3 M3 
NC_087233.1     ALN     nucleotide_to_protein_match     6521096 6521202 78      -       .       ID=match00005;Name=NC_087233.1_6523;Target=XP_051166906.1  16 52 +;Gap=M7 I2 M28 
NC_087233.1     ALN     nucleotide_to_protein_match     6520870 6521024 161     -       .       ID=match00005;Name=NC_087233.1_6523;Target=XP_051166906.1  53 104 +;Gap=M52 
NC_087233.1     ALN     nucleotide_to_protein_match     6520577 6520795 334     -       .       ID=match00005;Name=NC_087233.1_6523;Target=XP_051166906.1  105 177 +;Gap=M73 
NC_087233.1     ALN     nucleotide_to_protein_match     6520360 6520479 155     -       .       ID=match00005;Name=NC_087233.1_6523;Target=XP_051166906.1  178 217 +;Gap=M40 
NC_087233.1     ALN     nucleotide_to_protein_match     6520028 6520280 240     -       .       ID=match00005;Name=NC_087233.1_6523;Target=XP_051166906.1  218 301 +;Gap=M84 
NC_087233.1     ALN     nucleotide_to_protein_match     6519667 6519696 39      -       .       ID=match00005;Name=NC_087233.1_6523;Target=XP_051166906.1  302 311 +;Gap=M10 
NC_087233.1     ALN     nucleotide_to_protein_match     6519262 6519498 41      -       .       ID=match00005;Name=NC_087233.1_6523;Target=XP_051166906.1  312 394 +;Gap=M6 D1 M60 I4 M6 I1 M6 
NC_087233.1     ALN     nucleotide_to_protein_match     6519169 6519180 21      -       .       ID=match00005;Name=NC_087233.1_6523;Target=XP_051166906.1  395 399 +;Gap=M2 I1 M2 
NC_087233.1     ALN     nucleotide_to_protein_match     6518993 6519129 43      -       .       ID=match00005;Name=NC_087233.1_6523;Target=XP_051166906.1  400 446 +;Gap=M20 D1 M16 I1 M8 I1 M1 
NC_087233.1     ALN     nucleotide_to_protein_match     6518807 6518932 26      -       .       ID=match00005;Name=NC_087233.1_6523;Target=XP_051166906.1  447 486 +;Gap=M6 D1 M21 D1 M13 
NC_087233.1     ALN     nucleotide_to_protein_match     6518637 6518686 97      -       .       ID=match00005;Name=NC_087233.1_6523;Target=XP_051166906.1  487 503 +;Gap=M17 
NC_087233.1     ALN     nucleotide_to_protein_match     6518578 6518581 -40     -       .       ID=match00005;Name=NC_087233.1_6523;Target=XP_051166906.1  504 503 +;Gap=D1 

For example, please see this part where the end (2108443) is smaller than the start (2108444):

NC_087242.1     ALN     nucleotide_to_protein_match     2108444 2108443 5       +

And also, the 9th column of the same line (181 > 180):

181 180 +;Gap=

This pattern happens at the end of each alignment for these sequences.

Please let me know if you have an issue reproducing this result or anything else. Thanks again!

Cheers,
Dong-Ha

[ogotoh]

Dear Dong-Ha,

Such an anomaly can occur at the first or the last exon when a predicted codon does not align with any query residue. You may see what happens by running Spaln with -O1 option (alignment mode). For your example, the following command (-T option is added for potentially better performance)

$ spaln -Q7 -d gnm_test -T InsectAp -O1 ' toy_example.fa (2)'

will produces

>NC_087227.1 [1:256000]  ( 8151041 - 8407040 ) - >XP_023316908.1  [1:554]  ( 1 - 554 )
;C join(8278318..8278415,8278534..8278788,8278895..8279036,
;C 8279114..8279265,8279406..8279580,8279708..8279809,8279939..8280231,
;C 8280308..8280625,8280780..8280923,8282427..8282430)
PAM = 150, BIAS = 0.0, u = 2.0, v = 9.0
Score = 2365.8 (2380.7), 424.0 (=), 130.0 (#), 2.0 (g), 6.0 (u), (75.71 %)
ALIGNMENT   1 / 1

          M  T  L  P  F  H  G  T  E  P  R  K  K  E  E  I  L  V  H  A
 8278318 ATGACACTCCCCTTCCACGGGACTGAGCCGCGAAAAAAAGAGGAAATCCTCGTCCACGCG| NC_087227.1
       1  M  S  T  P  P  H  G  T  E  V  R  K  S  N  E  I  L  E  H  A | XP_023316908.1

          K  D  F  L  D  Q  Y  F  T  S  I  R  R
 8278378 AAGGACTTCCTGGACCAATATTTCACGTCCATTCGGAGgtgagttaatcattcagtgaat| NC_087227.1
      21  K  D  F  L  N  Q  Y  F  T  S  I  K  R                      | XP_023316908.1

...

          T  F  G  P  L  J  N  V  R
 8280898 ACATTCGGACCGTTGAGTAACGTTAGgttcgctgtatttgctctcggctcatcggcctat| NC_087227.1
     546  N  F  G  P  L  S  N  V  R                                  | XP_023316908.1

;; skip 1440 nt's

                                        O
 8282398 ggagagccactctacgtgtttgtccgcagGTAG                           | NC_087227.1
     555                                -                            | XP_023316908.1

Actually, I have no idea how to properly represent this situation in Gff3 format. In the default (-O4) output format, you can see that the genomic sequence does not match any part of the query.

I guess one possibility that the query amino acid sequence is not full length, lacking a C-terminus region, so that Spaln forcibly finds a nearby termination codon in this example.

I thought that -LS or -LC (local similarity) option can prevent such anomalies, but the current implementation does not work as expected. I will check the behavior of Spaln when -LS (or -LC) option is given.

Osamu,