conkit-plot peval fails to match sequences

[sadiogo]

General Information

• ConKit version: 0.12.0
• Python version: 3.6.9
• Environment (if applicable): WSL

Example

A minimal example to reproduce the error:

../.local/bin/conkit-plot cmap conkit.jones jones conkit.rr casprr
or
../.local/bin/conkit-plot cmap conkit.jones jones conkit.mat ccmpred

The .jones, .rr and .mat file were generated from the conkit-predict script, which worked fine. I used them to exclude the possibility that my original files were the problem; my original goal was use a reference structure to determine which contacts were true positives, but I was getting a traceback. So I tried just plotting without the reference structure and I got the same traceback.

Traceback

The Python traceback

File "/usr/lib/python3.6/runpy.py", line 193, in _run_module_as_main
    "__main__", mod_spec)
  File "/usr/lib/python3.6/runpy.py", line 85, in _run_code
    exec(code, run_globals)
  File "/home/sadiogo/.local/lib/python3.6/site-packages/conkit/command_line/conkit_plot.py", line 441, in <module>
    main()
  File "/home/sadiogo/.local/lib/python3.6/site-packages/conkit/command_line/conkit_plot.py", line 283, in main
    con.set_sequence_register()
  File "/home/sadiogo/.local/lib/python3.6/site-packages/conkit/core/contactmap.py", line 503, in set_sequence_register
    c.res1 = self.sequence.seq[res1_index - 1]
  File "/home/sadiogo/.local/lib/python3.6/site-packages/conkit/core/contact.py", line 268, in res1
    self._res1 = Contact._set_residue(amino_acid)
  File "/home/sadiogo/.local/lib/python3.6/site-packages/conkit/core/contact.py", line 483, in _set_residue
    raise ValueError("Unknown amino acid: {} (assert all is uppercase!)".format(amino_acid))
ValueError: Unknown amino acid: - (assert all is uppercase!)

[FilomenoSanchez]

Hi @sadiogo, I think the problem here might be that you are providing a MSA file (conkit.jones) instead of a FASTA file with the sequence of the protein of interest. Have you tried:

../.local/bin/conkit-plot cmap conkit.fasta fasta conkit.rr casprr

Where conkit.fasta is the FASTA file that you used as input for conkit-predict. Please let me know whether this solves the issue or the error persists.

[sadiogo]

Hi @FilomenoSanchez, I did try that. I get the following traceback when I use a single fasta sequence (ungapped) :

  File "/usr/lib/python3.6/runpy.py", line 193, in _run_module_as_main
    "__main__", mod_spec)
  File "/usr/lib/python3.6/runpy.py", line 85, in _run_code
    exec(code, run_globals)
  File "/home/sadiogo/.local/lib/python3.6/site-packages/conkit/command_line/conkit_plot.py", line 441, in <module>
    main()
  File "/home/sadiogo/.local/lib/python3.6/site-packages/conkit/command_line/conkit_plot.py", line 283, in main
    con.set_sequence_register()
  File "/home/sadiogo/.local/lib/python3.6/site-packages/conkit/core/contactmap.py", line 503, in set_sequence_register
    c.res1 = self.sequence.seq[res1_index - 1]
IndexError: string index out of range

When I use a single fasta sequence, but gapped, I get the same Unknown amino acid: - error.

By the way, I provided my own alignment when I used conkit-predict to generate the .mat and .rr files. Does this affect anything?

[FilomenoSanchez]

The traceback you have shown in your last comment occurs when the residue numbering in the contact map and the sequence do not match, particularly when the sequence in the FASTA file is shorter than it should. For example, the sequence in the FASTA file has 100 residues but the residue numbers in the contact map go all the way up to 200. I think this might have occurred because you provided a gapped alignment as input for conkit-predict, and the numbering in the resulting map was shifted because of these gaps. This might be something I'll have to fix later, but for the time being I think you should be able to fix this doing the following:

1. Take the single sequence FASTA file (but gapped) and replace gap characters "-" with "X", which stands for unknown residue.
2. Save the file conkit_modified.fasta and repeat the command:
   ../.local/bin/conkit-plot cmap conkit_modified.fasta fasta conkit.rr casprr

[sadiogo]

Thanks, that worked!

In retrospect, I had a gapped fasta alignment which I converted to a3m using conkit-convert, and that is what I used as input for conkit-predict. I opened the a3m file now and saw that the first sequence has gaps in it, which is not the way a3m files usually work, right? Perhaps there is a bug when converting fasta alignments to a3m?

[FilomenoSanchez]

I am not sure whether A3M format specifications explicitly prohibit gaps in the first sequence of the alignment? In any case, it is not trivial to remove gaps from the template sequence when converting from FASTA to A3M, as this would require to recalculate the alignment with all the sequences included in the MSA all over again (AFAIK), and it's well beyond conkit's original purpose. I think the only realistic fix here would be to make conkit interpret gaps in the FASTA sequence just as if they were unknown residues.

[sadiogo]

My knowledge on a3m files is entirely based on my experience using the MPI Bioinformatics Toolkit. From what I've seen, the first sequence is always gapless (which is propitious, since it is generally the query sequence). In any case, they have a service there called FormatSeq which can convert fasta alignments to a3m. I'm not sure the source code for that is available in the hhsuite, but you can reference people to the MPI Bioinformatics Toolkit if need be.

[FilomenoSanchez]

Yes indeed it seems that the A3M includes a gapless query sequence (my knowledge about this format is also limited). The hhstuite provides a script reformat.pl capable of dealing with this, it shouldn't be too difficult to reverse-engineer it into python and add it to conkit, I will open another issue for this. Since your original problem seems to be fixed now, I will close this issue but feel free to re-open it if otherwise.

[sadiogo]

One last question, @FilomenoSanchez. When plotting the reference structure to assess false positive contacts, the structure being added has to satisfy one or both of these conditions (?):

1. The number of residues in the PDB must be equal to the query sequence length?
2. The PDB file must be a structure model of the query sequence, i.e., they must have identical sequences?

I think these questions also apply to ConPlot? (very nice server, btw)

[FilomenoSanchez]

@sadiogo

In the case of ConPlot:

1. The number of residues does not need to be identical, what is important is that the sequence has at least as many residues as the PDB file.
2. The residue numbers in the PDB structure have to match the order of residues in the provided sequence. This way we account for gaps in the protein model. For example, for the following FASTA sequence:

 > sequence
MAVFY

You would need a PDB with the following residue numbers:

 1   2   3   4   5
MET ALA VAL PHE TYR

However, a PDB with missing residues 2-3 would still be valid, as long as the residue numbering is consistent:

 1   4   5
MET PHE TYR

Also please note that ConPlot cannot interpret gaps in the FASTA sequence. Conkit is a bit more clever dealing with gaps in PDB files, but I would generally recommend following the logic explained here when working with any of these two tools.

[sadiogo]

Okay, so from what you told me, this should work(?):

>sequence
MAVFYDD

PDB:

 1   4   5 
ILE TYR TRP

Which is equivalent to saying:

SEQ: MAVFYDD
PDB: I--YW--

[FilomenoSanchez]

Yes, this should work (if you have noticed it doesn't work in ConPlot, please open an issue at rigdenlab/conplot). Is the residue sequence mismatch intentional (residue 1 is MET in the FASTA and ILE in PDB)? The residue names are parsed from the FASTA sequence and are used in the tooltip display that appears when you hover over contacts and the diagonal of the plot (in ConPlot), so that information won't match the PDB (however this input will still produce a contact map).