wrong alphabet exception about jstacs HOT 6 CLOSED

Hi Rebecca,

thanks a lot for your email. I haven't seen this before.

The error clearly states that the amino acid "J" is unknown in our amino acid alphabet.
"J" encodes for leucine or isoleucine (https://www.ddbj.nig.ac.jp/ddbj/code-e.html).
It seems that "J" is used in a result of the search algorithm that you use (cf. "projects.gemoma.GeMoMa.addHit(GeMoMa.java:1173)"), which seems to be mmseqs. Could you please write which version you use?

Typically, broken GeMoMaPipeline runs leave some information in the directory "GeMoMa_temp". Here you could check the results of the search algorithm identifying the gene model that has matches including the amino acid "J". This could help to create a minimal example.

GeMoMa uses the amino acids given in the substitution matrix to define the amino acid alphabet. The substitution matrix that we currently use do not have entries for "J" (https://github.com/Jstacs/Jstacs/blob/master/projects/gemoma/test_data/BLOSUM62.txt). This matrix was downloaded from NCBI some time ago to be compatible to blast (ftp://ftp.ncbi.nlm.nih.gov/blast/matrices/BLOSUM62).

Using a substitution matrix that has entries for "J" might solve the problem. I just had a look and found: https://www.ncbi.nlm.nih.gov/IEB/ToolBox/C_DOC/lxr/source/data/BLOSUM62
Based on this, we might be able to create such a substitution matrix in the correct format.

However, I compared our BLOSUM matrix with the one above and found that they differ not only in the column for "J" but also for some individual values. Hence, it would be good to know which substitution matrix was used by the search algorithm.

Best regards, Jens

from jstacs.

Thanks so much for the reply! I am trying to dig through the temp files. I uploaded them here (along with a log file) in case you have time to have a look: https://www.dropbox.com/s/w9a8mf21vvurbdz/gemoma-test.tar.gz?dl=0

Thanks again!

from jstacs.

Hi Rebecca,

thanks a lot for the huge amount of data.

I checked the log file and found that you are using mmseqs version 13.45111.
Checking mmseqs substitution matrices: https://github.com/soedinglab/MMseqs2/tree/master/data
I did not find a column for amino acid "J" in BLOSUM62 and some others that I checked.
I also did not find that they have changed the default substitution matrix or anything regarding leucine/isoleucine: https://github.com/soedinglab/MMseqs2/releases

In addition, I had a look into GeMoMaPipeline-16045397510071019339/0/mmseqs.tabular. Only 3 out of more then 12M lines contained a "J" in the alignment:

3573076 gene_16314_1 NW_021940722.1 44.400 27 15 0 15 41 999 1079 2.172E+01 29 - 62 12 - - - PELPSGAGPFDAATTPARLSDATAFPP PRVPSGPSPRSVSLJPDRRCQATALPP 61 4478
3850325 gene_16314_3 NW_021940722.1 44.400 27 15 0 14 40 999 1079 2.078E+01 29 - 62 12 - - - PELPSGAGPFDAATTPARLSDATAFPP PRVPSGPSPRSVSLJPDRRCQATALPP 60 4478
8552716 gene_16474_0 NW_021939122.1 35.800 67 39 3 11 75 11486396 11486202 8.653E+00 32 - 70 24 - - - LLARLGLQALLALLPALARGLGLHPRPPCPGPGSAAPALPAGRAAPLWGPAR--GDEYERRYSRAFP MLRAVGLNAGRG-APPJSRH-GTLPFPGAPGPARAAPECPSPRSGTERGPISPVGEQPESRENKYGP 121 17883582

(first column is the line number starting with 0)

A very simple work around is to remove the corresponding reference transcripts to avoid "J" in the search results of mmseqs. However, I would recommend to contact the mmseqs developers. Maybe they can give better advices.

Looking forward to any news.

best regards, Jens

from jstacs.

I've implemented a new paremeter "replace unknown" that allows to replace unknown amino acids in the search results by "X".

c93cc06

Will be available with the next release.

from jstacs.

Thank you so much!

from jstacs.

Laura Baxter also reported this problem to me. However, she found another way to avoid this problem:

"I managed to track the source of the issue to nucleotide ambiguity codes in the target genome sequence FASTA file. I replaced all non-ATGC characters in the sequence with N with a perl one-liner, and then everything worked great. Worth remembering in case anyone else has the same issue in the meantime.

The one-liner to replace any non-ATGC character in the sequence with N was:
perl -pe 's/[^AGTC\n]/N/gi unless m/>/;' oldfile.fas > newfile.fas
"

Thanks Laura.

from jstacs.